Optimization with trapz and data

The basic bit. I have an experimental data point. I can estimate what the value should be from physics; the value is dependent on temperature(a constant) and energy(covers a range). I want to find the value of temperature that (chi) minimizes the experimental value to my estimated value to find the temperature.

The complex bit: I have to integrate over energy to get my estimated signal. The integral is non-trivial. It involves a predined function (chosen based on physics). Pretty much this is an exponential with energy and temperature outfront and in the exponent.

This is multipled by two emperical functions that take into account filters and diagnostic response functions that depend on energy. The emperical functions aren't really functions; I have emperical data, did multiple fits to sections of the data to get the best fit possible and compiled it all into an array using the relevant energy values. A single fit is not possible.

Currently I'm solving this with for loops, generating large amounts of data. There are actually 9 data points, which should be solved consistently with each other. It takes 10-20 minutes with a parallel loops to run the current iteration. The issue is, we think there's a second temperature. Running it with loops was already inefficient but I've never used the optimziation toolbox and kept getting stuck trying to use it. A second temperature makes doing this with loops....pretty much impossible.

I'm wondering if there's a way to use the optimization toolbox (or something else in Matlab?) to improve solving this problem.

[darova]

Can you the script? Some images?

[plasmageek]

Hmm...most of it won't be useful I think. This is the main workhorse of the process

%%Assuming this functional form for distribution

Max2F = 2*sqrt(E)*(exp(-E/kT1)./(pi*kT1)^(3/2)+exp(-E/kT2)./(pi*kT2)^(3/2));

%%Generate 'function' to 'integrate' ; take area under curve

intFnc = Max2F.*Trans.*IP_response'*Omega/A;

EstVal = trapz(E,intFnc);

%%Compare estimated value to exerpimental value

Chi = (EstVal-ExpVal)^2/ExpVal;

I'm trying to find the temperatures that minimize my Chi function. Trans and IP_Response are emperical data. Omega and A are constants.

[darova]

I think fit() will be good

[plasmageek]

Hmm I'm not sure I follow. Normally I have a complete set of x-y data to fit, but I don't really have that in this case.

[plasmageek]

So here's my current issue. I need to solve all my channels consistently. I was able to combine my emperical data sets and find a fit for them. The issue is this fit changes for each channel, but it has to be included in the function. I'm trying to use lsqcurvefit. Below is my current code (which doesn't run). My xdata represents my different channels.

The variable E is an array from 1 to 300. Omega is a single value for each channel.

Area   = a*1e-3*b*1e-3; %%m^2
for ii = 1:7
    idx    = ii+1;
    const  = T(idx,:).*IP_response';
    
    [E,const] = prepareCurveData(E,const);
    Ft{ii}    = fit(E,const,'smoothingspline');
end
%%
xdata = 2:8; %%calling out channels of interest
ydata = double(SigFin(xdata))/IP_Fade;

fun = @(x,xdata)Omega(xdata)/Area*trapz((2*sqrt(x(1))/(x(2)^(3/2)))*exp(-E/x(2)).*Ft{xdata}(E));
x0  = [1, .020];
X   = lsqcurvefit(fun,x0,xdata,ydata);

[plasmageek]

And my other thought was to make a function and have it called. I put together the following, but it's obviously not right either...

function fun = HerieAnalysis(x,Area,Omega,xdata,E,Ft)

for ii = 1:length(xdata)
    switch xdata(ii)
        case 1
            const = Ft{1};
            Ang   = Omega(1);
        case 2
            const = Ft{2};
            Ang   = Omega(2);
        case 3
            const = Ft{3};
            Ang   = Omega(3);
        case 4
            const = Ft{4};
            Ang   = Omega(4);
        case 5
            const = Ft{5};
            Ang   = Omega(5);
        case 6
            const = Ft{6};
            Ang   = Omega(6);
        case 7
            const = Ft{7};
            Ang   = Omega(7);
        case 8
            const = Ft{8};
            Ang   = Omega(8);
        case 9
            const = Ft{9};
            Ang   = Omega(9);
    end
    fun = Ang/Area*trapz((2*sqrt(x(1))/(x(2)^(3/2)))*exp(-E/x(2)).*const(E));
end