2013-05-25 6 views
5

मैं अपने मैटलैब का उपयोग कर रहा हूं, लेकिन अंततः पाइथन में मेरे सभी विश्लेषण करने के लिए यह मेरी दृष्टि है क्योंकि यह एक वास्तविक प्रोग्रामिंग भाषा है और कुछ अन्य कारण हैं।कम स्क्वायर न्यूनतमकरण कॉम्प्लेक्स नंबर

हाल ही की समस्या जो मैं निपटाने का प्रयास कर रहा हूं वह जटिल डेटा के कम से कम वर्गों को कम करने के लिए है। मैं एक इंजीनियर हूं और हम अक्सर जटिल प्रतिबाधा से निपटते हैं, और मैं मापा डेटा के लिए एक सरल सर्किट मॉडल फिट करने के लिए वक्र फिटिंग का उपयोग करने की कोशिश कर रहा हूं।

प्रतिबाधा समीकरण इस प्रकार है:

जेड (डब्ल्यू) = 1/(1/आर + j * डब्ल्यू * सी) + j * डब्ल्यू * एल

मैं तो खोजने की कोशिश कर रहा हूँ आर, सी, और एल के मूल्य जैसे कि कम से कम वर्ग वक्र पाया जाता है।

मैंने ऑप्टिमाइज़ेशन पैकेज जैसे ऑप्टिमाइज़ेशन.curve_fit या optimize.leastsq का उपयोग करने का प्रयास किया है, लेकिन वे जटिल संख्याओं के साथ काम नहीं करते हैं।

तब मैंने अपने अवशिष्ट कार्य को जटिल डेटा की परिमाण को वापस करने की कोशिश की, लेकिन यह या तो काम नहीं किया।

+2

यह एक आशाजनक उत्तर के साथ बिल्कुल वही समस्या प्रतीत होता है: http://stackoverflow.com/questions/14296790/python-scipy-leastsq-fit-with-complex-numbers – jmihalicza

उत्तर

4

अंत में, लक्ष्य मॉडल और के बीच चुकता मतभेद की राशि का निरपेक्ष मान मनाया Z को कम करना है:

abs(((model(w, *params) - Z)**2).sum()) 

मेरे original answer एक residuals समारोह जो एक अदिश रिटर्न के लिए leastsq लागू करने का सुझाव दिया वास्तविक और काल्पनिक मतभेदों के वर्गों का योग दर्शाता है:

def residuals(params, w, Z): 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    return diff.real**2 + diff.imag**2 

Mike Sulzer suggested अवशिष्ट फ़ंक्शन का उपयोग करके जो फ़्लोट्स का वेक्टर लौटाता है।

from __future__ import print_function 
import random 
import numpy as np 
import scipy.optimize as optimize 
j = 1j 

def model1(w, R, C, L): 
    Z = 1.0/(1.0/R + j*w*C) + j*w*L 
    return Z 

def model2(w, R, C, L): 
    Z = 1.0/(1.0/R + j*w*C) + j*w*L 
    # make Z non-contiguous and of a different complex dtype 
    Z = np.repeat(Z, 2) 
    Z = Z[::2] 
    Z = Z.astype(np.complex64) 
    return Z 

def make_data(R, C, L): 
    N = 10000 
    w = np.linspace(0.1, 2, N) 
    Z = model(w, R, C, L) + 0.1*(np.random.random(N) + j*np.random.random(N)) 
    return w, Z 

def residuals(params, w, Z): 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    return diff.real**2 + diff.imag**2 

def MS_residuals(params, w, Z): 
    """ 
    https://stackoverflow.com/a/20104454/190597 (Mike Sulzer) 
    """ 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    z1d = np.zeros(Z.size*2, dtype=np.float64) 
    z1d[0:z1d.size:2] = diff.real 
    z1d[1:z1d.size:2] = diff.imag 
    return z1d 

def alt_residuals(params, w, Z): 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    return diff.astype(np.complex128).view(np.float64) 

def compare(*funcs): 
    fmt = '{:15} | {:37} | {:17} | {:6}' 
    header = fmt.format('name', 'params', 'sum(residuals**2)', 'ncalls') 
    print('{}\n{}'.format(header, '-'*len(header))) 
    fmt = '{:15} | {:37} | {:17.2f} | {:6}' 
    for resfunc in funcs: 
     # params, cov = optimize.leastsq(resfunc, p_guess, args=(w, Z)) 
     params, cov, infodict, mesg, ier = optimize.leastsq(
      resfunc, p_guess, args=(w, Z), 
      full_output=True) 
     ssr = abs(((model(w, *params) - Z)**2).sum()) 
     print(fmt.format(resfunc.__name__, params, ssr, infodict['nfev'])) 
    print(end='\n') 

R, C, L = 3, 2, 4 
p_guess = 1, 1, 1 
seed = 2013 

model = model1 
np.random.seed(seed) 
w, Z = make_data(R, C, L) 
assert np.allclose(model1(w, R, C, L), model2(w, R, C, L)) 

print('Using model1') 
compare(residuals, MS_residuals, alt_residuals) 

model = model2 
print('Using model2') 
compare(residuals, MS_residuals, alt_residuals) 

पैदावार

Using model1 
name   | params        | sum(residuals**2) | ncalls 
------------------------------------------------------------------------------------ 
residuals  | [ 2.86950167 1.94245378 4.04362841] |    9.41 |  89 
MS_residuals | [ 2.85311972 1.94525477 4.04363883] |    9.26 |  29 
alt_residuals | [ 2.85311972 1.94525477 4.04363883] |    9.26 |  29 

Using model2 
name   | params        | sum(residuals**2) | ncalls 
------------------------------------------------------------------------------------ 
residuals  | [ 2.86590332 1.9326829 4.0450271 ] |    7.81 | 483 
MS_residuals | [ 2.85422448 1.94853383 4.04333851] |    9.78 | 754 
alt_residuals | [ 2.85422448 1.94853383 4.04333851] |    9.78 | 754 

तो ऐसा लगता है मॉडल समारोह पर निर्भर हो सकता उपयोग करने के लिए जो अवशिष्ट समारोह:

यहाँ इन अवशिष्ट कार्यों का उपयोग कर परिणाम की तुलना है। model1 और model2 की समानता के परिणामस्वरूप अंतर में व्याख्या करने के लिए मैं हूं।

+0

जटिल संख्या विभिन्न चौड़ाई वाली फ्लोट पर आधारित हो सकती है, numpy दस्तावेज़ीकरण के आधार पर यह दृश्य डेटा के बाइट लेवल पुनर्वितरण को निष्पादित करता है, यानी फ्लोट np.float32 (reres भी float16 और float64 को refereeing को समाप्त कर सकता है, np.complex128 का उपयोग करते समय, इसे 4 फ्लोट में परिवर्तित किया जाएगा, और वास्तविक फ्लोटिंग पॉइंट नंबरों की संरचना के कारण मानों का अर्थ कम होगा, माइक Sulzer का जवाब सुरक्षित है, क्योंकि मान सही फ़्लोटिंग पॉइंट आकार में परिवर्तित हो जाएगा। आप भी निरंतर स्मृति, यानी आइटम आकार = 16 बाइट्स, लेकिन strides = 20 बाइट्स, यह अभी भी एक 1 डी सरणी है। – glenflet

+0

@glenflet: सुधार के लिए धन्यवाद। dtype mismatch समस्याद्वारा तय की जा सकती है 210 बदल रहा है 'diff.view (' float ') ' ' diff.astype (np.complex128) .view (np.float64) '। (मैंने ऊपर दिए गए उत्तर को पर संपादित किया है।) गैर-संगतता समस्या उत्पन्न नहीं होती है क्योंकि 'diff' दो सरणी का अंतर है:' diff = model (w, आर, सी, एल) - Z' । NumPy हमेशा 'मॉडल (डब्ल्यू, आर, सी, एल) - जेड और के परिणाम वाले एक संगत सरणी बनायेगा, पाइथन वैरिएबल नाम' diff' को बांधता है। – unutbu

+1

मैं कहूंगा कि आपके पास दो मॉडलों के बीच का अंतर है फ्लोटिंग पॉइंट 'अवशेष' का आकार मॉडल 2 के लिए फ्लोट 32 का उपयोग करता है क्योंकि जेड np.complex64 यानी 2 32 बिट फ्लोट्स है, अन्य दो फ़ंक्शन फ़्लैट 64 पर कास्टिंग कर रहे हैं, आपको तुलना करनी चाहिए फ़ोकस कॉल की संख्या भी, डॉक्टर के मान से मेरा मानना ​​है कि स्टेप साइज 2 * sqrt (मशीन परिशुद्धता) के लिए डिफ़ॉल्ट है, जो फ्लोट 32 के लिए बड़ा है, यदि आप maxfev (आपके कोड में 800) से अधिक है तो यह महत्वपूर्ण है। इसके अलावा यह एक स्थानीय मिम पर छोड़ सकता है, या वर्गों के योग में उनकी त्रुटि हो सकती है, हालांकि इसका कोई प्रभाव नहीं होना चाहिए। – glenflet

5

to unutbu answer's का जिक्र करते हुए, कार्य अवशिष्ट में परिमाण वर्ग को प्राप्त करके उपलब्ध जानकारी को कम करने की आवश्यकता नहीं है क्योंकि leastsq परवाह नहीं है कि संख्या वास्तविक या जटिल हैं, लेकिन केवल वे ही 1 डी सरणी के रूप में व्यक्त की जाती हैं कार्यात्मक संबंध की अखंडता।

def residuals(params, w, Z): 
    R, C, L = params 
    diff = model(w, R, C, L) - Z 
    z1d = np.zeros(Z.size*2, dtype = np.float64) 
    z1d[0:z1d.size:2] = diff.real 
    z1d[1:z1d.size:2] = diff.imag 
    return z1d 

यह केवल परिवर्तन किए जाने की जरूरत है:

यहाँ एक प्रतिस्थापन बच कार्य है। बीज 2013 के उत्तर हैं: [2.96564781, 1.99929516, 4.00106534]।

to unutbu answer's से संबंधित त्रुटियों को परिमाण के क्रम से काफी अधिक कम किया जाता है।

संबंधित मुद्दे