Like Tree9Likes

Parameter estimation Nelder Mead

This is a discussion on Parameter estimation Nelder Mead within the C Programming forums, part of the General Programming Boards category; I am trying to minimize a function for three parameters. The code for main and the function I am trying ...

  1. #1
    Registered User
    Join Date
    Apr 2008
    Posts
    190

    Parameter estimation Nelder Mead

    I am trying to minimize a function for three parameters. The code for main and the function I am trying to minimize are below. It runs but doesnt get to the correct solution. I am running the same thing in Matlab and it gets it right.

    Just wondering if anyone can see a reason it may not converge as accurately as the equivelent matlab.

    Thanks

    Alex

    Code:
    double Actual_Output[5];
    double Input[5];
    double a_global, b_global;
    //double x_coeff_global;
    float init_T[23][2]={0};
    float eta_glob, phi_glob, omega_2_glob, c1_glob;
    float Root(float left, float right, float tol, int *count);
    int num;
    float f(float x);
    //Function to find physical parameters a and b (stored globally as a_global and b_global)
    
    main()
    {
        float left, right, tol, fall_angle, fall_time, guess, c0, xx;
    	int	n, i,num_inputs = 23; //number of parameters to search for
    	double	x[3], y[3];
    	double	length = 1.00;
    	double	fopt;
    	int	timeout = 10000;
    	double	eps = 1.0e-12;//accuracy of Nelder Mead
    	left = 0;
        tol = 0.0001;//bisection accuracy
        FILE *fp;
    //Starting estimates for parameters a and b (Check documentation (EQN 40))
    	x[0] = 0.5;//a
    	x[1] = 1;//b
    	x[2] = -1;//b
    	//Timing Values for 5 revolutions
        Actual_Output[0] = 1.2;
        Actual_Output[1] = 2.693;
        Actual_Output[2] = 4.325;
        Actual_Output[3] = 6.131;
        Actual_Output[4] = 8.125;
    
        //Revolution Number corresponding with revolution time above
        Input[0] = 1;
        Input[1] = 2;
        Input[2] = 3;
        Input[3] = 4;
        Input[4] = 5;
        n=3;
        //calculating physical condition parameters a and b (stored as a_global and b_global)
    	if (NelderMeadSimplexMethod(n, function, x, length, &fopt, timeout, eps) == success) {
    		printf("reaching to minimum ");
    	} else {
    		printf("timeout  ");
    	}
    	printf(" [ %lf %lf %lf]  %lf\n", x[0], x[1],x[2], fopt);
    	//At this point a_global and b_global are assigned. These are required for fall parameters to calculate c1:
    	//beta:
    	//beta = a_global*b_global*b_global; 
    return 0;
    }
    
    static double function(int n, double x[])
    {
    	double c;
        double Fitted_Curve[5];
        double Error_Vector[5];
        int i;
        double a, b, sum = 0;
        
    //    Actual_Output[5] = {1.2, 2.693, 4.325, 6.131, 8.125};
    
    a = x[0];
    b=x[1];
    c=x[2];
    
      for (i = 0; i <= 4; i++)
      {
    
        Fitted_Curve[i] = (1/(a*b))*(c-asinh(sinh(c)*exp(a*Input[i]*2*pi)));
    
        Error_Vector[i] = Actual_Output[i]-Fitted_Curve[i];
        }
    
        for (i = 0; i <= 4; i++)
        {
            sum = sum + Error_Vector[i]*Error_Vector[i];
        }
        printf("sum = %f\n", sum);
        a_global = a;
        b_global = b;
       // x_coeff_global = x_coeff;
        return sum;
    
    }
    Last edited by a.mlw.walker; 08-16-2011 at 11:04 AM.

  2. #2
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    I'm going to give you the same advice here as before: do it yourself with a calculator, and print it out all the intermediate results -- since it seems that the very first iteration is wrong(?), then you don't even have to loop; just copy the code into a new file, remove the loop to just do one iteration and look at it. Find the line that doesn't calculate the way it should, and fix it.
    AndrewHunter likes this.

  3. #3
    Registered User
    Join Date
    Apr 2008
    Posts
    190
    Ok, that is the solution isnt it. If it works its not a C problem so I will just do it by maths...

  4. #4
    Registered User
    Join Date
    Nov 2010
    Location
    Long Beach, CA
    Posts
    5,443
    You need to tell us exactly what's incorrect. "Wrong answer" doesn't help us much. Are you off by .01%? Are you off by 1000%? Also, as tabstop pointed out, determine which step is wrong by printing out after every step in the C program (or examining variables while stepping through the code in a debugger) and every step in the Matlab program (this is default behavior, unless you append a ; to your matlab statements) and comparing. Then you can isolate the problem to a small calculation, and if you still can't find it, we will have a much easier time helping you.

    Another thing that was brought up before is to use doubles instead of float. Matlab uses doubles by default, so you should always use doubles in your C code to avoid any rounding/truncation/precision errors. Make sure you always use the double version of functions (fabs, etc).

    Also, a proper C program is:
    Code:
    int main(void)
    {
        return 0;
    }
    Always be explicit in what parameters you expect and what you return, especially for main. Read why it matters here: Cprogramming.com FAQ > main() / void main() / int main() / int main(void) / int main(int argc, char *argv[])

    Lastly, your use of global variables worries me. There's always the potential for shadowing, difficulty of determining scope and difficulty of debugging when you use globals. Read this brief article and see if it convinces you to move those things into functions and pass them around. You can also group related variables in a struct and pass that around instead of the individual parts to help keep things terse.

  5. #5
    Registered User
    Join Date
    Apr 2008
    Posts
    190
    OK so I have used my calculator to get the results that it is calculating for the first iteration of ht efunction. a = 1, b = 1, c = -1. The image attached shows the results for the below code. It is exactly what my calculator gets and also the matlab (both attached). I just think it is stopping early or soemthing...
    the console is the c for the first iteration of the function, the other image is the matlab output - the sse value is just the first one, not the total sum...
    Attached Images Attached Images   

  6. #6
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    Quote Originally Posted by a.mlw.walker View Post
    OK so I have used my calculator to get the results that it is calculating for the first iteration of ht efunction. a = 1, b = 1, c = -1.
    Are you sure that's what you want? Your code says
    Code:
    x[0] = 0.5;//a
    I hate to ask, but: you do realize you never actually call the function in your code? Check what NelderMeadSimplexMethod does, because who knows what it does.

  7. #7
    Registered User
    Join Date
    Apr 2008
    Posts
    190
    The matlab starts with the initial values i posted, thats why i used them. that 0.5 is closer to the answer, i changed it to that to see if my guess was too far away. I have checked Nelder Mead it calls the function and that is how i produced those values in the console. i didnt remove it from the code, just put a break point to stop it computing any further. I think it is ending the iterations early of something but the tolerance doesnt seem to affect the results really.

  8. #8
    Registered User
    Join Date
    Apr 2008
    Posts
    190
    Sorry anduril462, didnt see your post. I have checked the maths as you can see in my last post it is computing the same thing as the matlab, however the final results is different. a = 1, b = 1, c = -1 are the starting values for both, matlab returns as final results:
    a = 0.0250 b= 1.8944 c= -0.3690
    C code produces
    a = 3.52, b = 3.634 c= -0.000073
    i.e differrent

  9. #9
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    Quote Originally Posted by a.mlw.walker View Post
    The matlab starts with the initial values i posted, thats why i used them. that 0.5 is closer to the answer, i changed it to that to see if my guess was too far away. I have checked Nelder Mead it calls the function and that is how i produced those values in the console. i didnt remove it from the code, just put a break point to stop it computing any further. I think it is ending the iterations early of something but the tolerance doesnt seem to affect the results really.
    Did you inadvertently put the breakpoint somewhere in your loop of adding up the sum (that seems to be the only way that you get that answer out of your code)?

    Also, remember the warning about global variables: you can't assume that any of those global variables are the same from when you call NelderMead to when NelderMead calls your function, because NelderMead (or any function that it calls before it calls your function) can change any or all of those global variables at will.

  10. #10
    Registered User
    Join Date
    Apr 2008
    Posts
    190
    In terms of the global variables, I have tried to be careful that I dont have that problem however I understand the problem. I never managed to pass variables suvvessfully through the NelderMead function, so I just went global.
    Sum is correct for the first iteration

  11. #11
    Registered User
    Join Date
    Apr 2008
    Posts
    190
    Ok I have attached three iterations from the c and matlab, printing everything...
    Attached Images Attached Images   

  12. #12
    ATH0 quzah's Avatar
    Join Date
    Oct 2001
    Posts
    14,826
    Posting screen shots of numbers doesn't really do anything for us here. As a whole, we don't really care about your math. We care about the C aspect of your program. If it has correct syntax, and doesn't invoke udefined behavior, then it's up to you to sort out your math issue.

    Honestly, I doubt there are many people around here that want to reverse engineer your set of numbers to figure out what math you have done to them.


    Quzah.
    Hope is the first step on the road to disappointment.

  13. #13
    Registered User
    Join Date
    Apr 2008
    Posts
    190
    This is a c issue. I have shown the mathematical side is correct in an earlier post. But C is doing something meaning that future calculations are not being computed correctly - I'm not sure what though - its a C thing not a math thing...

  14. #14
    ATH0 quzah's Avatar
    Join Date
    Oct 2001
    Posts
    14,826
    Quote Originally Posted by a.mlw.walker View Post
    its a C thing not a math thing...
    Prove it. You aren't showing your entire code.


    Quzah.
    Hope is the first step on the road to disappointment.

  15. #15
    Banned
    Join Date
    Aug 2010
    Location
    Ontario Canada
    Posts
    9,547
    Quote Originally Posted by a.mlw.walker View Post
    This is a c issue. I have shown the mathematical side is correct in an earlier post. But C is doing something meaning that future calculations are not being computed correctly - I'm not sure what though - its a C thing not a math thing...
    No offense is intended... but are you sure it's not an "I don't understand math in C" thing...

    If you're like most people you are trying to continue a previous experience into a new one... I did it when moving from Pascal to C wrote some horribly convoluted code in the process... In this case I get the feeling you're trying to make C behave like Matlab... which it can't do.
    quzah likes this.

Page 1 of 7 1234567 LastLast
Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Problem with simple Factorial Estimation
    By Arthurdent in forum C Programming
    Replies: 7
    Last Post: 06-16-2011, 07:28 PM
  2. Probability estimation
    By Mario F. in forum General Discussions
    Replies: 3
    Last Post: 09-18-2009, 08:40 AM
  3. Replies: 6
    Last Post: 01-08-2008, 09:25 AM
  4. project help (mead,median,mode)
    By Pliskin_Vamp in forum C Programming
    Replies: 1
    Last Post: 04-09-2007, 05:08 PM
  5. area estimation of graph
    By hei in forum C Programming
    Replies: 3
    Last Post: 10-16-2001, 12:23 AM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21