0

I am trying to get an average value for parameters to then plot with a given function. I think I have to somehow fill a 3-column array and then take the average of values of that array. I want to create 1000 values for popt[0] , popt[1] , and popt[2] and then take the average of all those values and then plot them.

for n in range(0,1000):
    params=np.zeros(3,1000)
    y3=y2+np.random.normal(loc=0.0,scale=0.1*y2)
    popt,pcov=optimize.curve_fit(fluxmeasureMW,bands,y3)
    params.append(popt[0],popt[1],popt[2])
    a_avg=st.mean(params[0:])
    b_avg=st.mean(params[1:])
    e_avg=st.mean(params[2:])

The final goal is to plot:

fluxmeasureMW(bands,a_avg,b_avg,e_avg)

I am just not sure how to iterate the fitting function to then output 1000 values. 1000 is arbitrary, I just want a good sample size. The values for y2 and bands are already defined and can be plotted without issue, as well as the function fluxmeasureMW.

ski
  • 7
  • 2
  • It isn't very clear what you are trying to do. Does `fluxmeasureMW` return a scalar or an array? Is `bands` defined - what is it? What should `a(b)(e)_avg` be - scalars or arrays? assuming `scipy.optimize.curve_fit`, `popt` are the optimized (non-dependent variable) arguments - why would you take their mean? – wwii Apr 11 '20 at 18:41
  • fluxmeasureMW returns an array, bands is an array of values that will never change, a_avg, b_avg, and e_avg should all be scalars. I want to take their mean because I have written the code such that every time it is run adds some noise to the data, but this noise changes every time the code is run so I keep getting different values for the parameters. – ski Apr 11 '20 at 18:52
  • You want to run the fit a thousand times adding a bit of noise to the `y` argument each time then get the mean of the optimized parameters.? There are three parameters being optimized? – wwii Apr 11 '20 at 23:46
  • Are `a(b)(e)_avg` positional or keyword arguments? `y2` is a scalar? – wwii Apr 11 '20 at 23:52

1 Answers1

0

Say your function is like this

def fluxmeasureMW(x,f,g,h):
    return result_of_calc

Just run the fit in a loop; accumulate the popts in a list then take the mean

from scipy import optimize
import numpy as np

n = 1000
t = []
for i in range(n):
    y3 = y2 + np.random.normal(loc=0.0,scale=.1*y2)
    popt,pcov = optimize.curve_fit(fluxmeasureMW,bands,y3)
    t.append(popt)

f,g,h = np.mean(t,0)

t will be a list of lists...

[[f,g,h],
 [f,g,h],
 ...]

np.mean(t,0) will average the values over the columns.

You could also use

import statistics
a = [[0, 1, 2],
     [1, 2, 3],
     [2, 3, 4],
     [3, 4, 5]]

for column in zip(*a):
    #print(column)
    print(statistics.mean(column))
wwii
  • 23,232
  • 7
  • 37
  • 77
  • Thank you! I will see if this works. Could you also please explain what np.mean(t,0) is doing? I understand it's taking the mean but I don't know why the zero is needed – ski Apr 12 '20 at 16:23
  • That did exactly what it was supposed to do, thank you! – ski Apr 12 '20 at 16:29
  • How do I go about plotting all the values in one column vs the values in another? Right now I have: `plt.plot(t[0],t[1],'o',label='a vs b')` `plt.plot(t[0],t[2],'o',label='a vs e')` `plt.plot(t[1],t[3],'o',label='b vs e')` – ski Apr 12 '20 at 18:09
  • `a,b,e = zip(*t)` will give you the three separate *columns* - it [transposes the list of lists](https://stackoverflow.com/questions/6473679/transpose-list-of-lists). Then `plt.plot(a,b);plt.plot(a,e);plt.plot(b,e)` if that was your intent. [Pyplot tutorial](https://matplotlib.org/tutorials/introductory/pyplot.html#sphx-glr-tutorials-introductory-pyplot-py) – wwii Apr 12 '20 at 18:51
  • Is there a way to make an average covariance matrix? I want to print the covariance matrix for the fitted parameters, but right now I can only print the last covariance matrix for the last fit in the for-loop. – ski Apr 20 '20 at 17:12
  • Searching with `average covariance matrix` produces results - some Q&A's in stats.stackexchange.com and stackoverflow.com – wwii Apr 20 '20 at 17:48