Monte Carlo Markov Chain - Metropolis-Hastings - Estimation of parameters











up vote
0
down vote

favorite












I have 6 parameters to estimate : $p=(theta=[a,b]$, $nu=[r_0,c_0,alpha,beta])$ with Bayesian and MCMC methods :



$$text{PSF}(r,c) = bigg(1 + dfrac{r^2 + c^2}{alpha^2}bigg)^{-beta}$$



and the modeling :



$$d(r,c) = a cdot text {PSF}_{alpha,beta} (r-r_0,c-c_0) + b + epsilon (r, c)$$



with $epsilon$ being a white gaussian noise.



This is the matrix form :



$$begin{bmatrix}d(1,1) \ d(1,2) \ d(1,3) \ vdots \ d(20,20) end{bmatrix}
= begin{bmatrix} text{PSF}_{alpha,beta}(1-r_0,1-c_0) & 1 \ text{PSF}_{alpha,beta}(1-r_0,2-c_0) & 1 \ text{PSF}_{alpha,beta}(1-r_0,3-c_0) & 1 \ vdots & vdots \ text{PSF}_{alpha,beta}(20-r_0,20-c_0) & 1 end{bmatrix} times begin{bmatrix}a \ b end{bmatrix}
+ begin{bmatrix}epsilon(1,1) \ epsilon(1,2) \ epsilon(1,3) \ vdots \ epsilon(20,20) end{bmatrix}$$



So we can write for the 1D data vector "d" :



$$d=H(nu),theta + epsilon$$ with $H$ the matrix defined above.



I know that we have the relation for posterior function ($d$=data and $p$ = vector of parameters :



$$f(p|d) = dfrac{f(p),f(d|p)}{int_{p}f(d|p),text{d}p}$$



We can also write :



$$f(p|d) propto text{Likelihood(d|p)},f(p)quad(1)$$ with $f(p)$ the prior function (that I can take as uniform distribution).



Now, how can I estimate the parameters $p$ from this relation $(1)$ with MCMC methods, especially in my case with Metropolis algorithm ?



With likelihood method, I used previously the following cost function to estimate these 6 parameters :



function cost = Crit_J(p,D)



% Compute the model corresponding to parameters p                                        

[R,C] = size(D);

[Cols,Rows] = meshgrid(1:C,1:R);                                                         

% Model

Model = (1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6));                             

model = Model(:);                                                                        

d = D(:);

% Introduce H matrix 

H = [ model, ones(length(model),1)];                                                     

% Compute the cost function : taking absolute value                                      

cost = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));                                      



end


And after, I perform these estimations with "Matlab fminsearch" function to find a local minimum.



Here, at first sight, I thought that I have to compute the distribution for each $p$ parameters generated randomly, but it seems to be more subtl.



Question 1) How can we proove that, starting from arbitrary values for parameters $p$, the Metropolis algorithm will converge to right estimations ?



Question 2) How to implement it for this concrete example ?



Any help is welcome, regards






maximum-likelihood monte-carlo log-likelihood







share|cite|improve this question


























    up vote
    0
    down vote

    favorite












    I have 6 parameters to estimate : $p=(theta=[a,b]$, $nu=[r_0,c_0,alpha,beta])$ with Bayesian and MCMC methods :



    $$text{PSF}(r,c) = bigg(1 + dfrac{r^2 + c^2}{alpha^2}bigg)^{-beta}$$



    and the modeling :



    $$d(r,c) = a cdot text {PSF}_{alpha,beta} (r-r_0,c-c_0) + b + epsilon (r, c)$$



    with $epsilon$ being a white gaussian noise.



    This is the matrix form :



    $$begin{bmatrix}d(1,1) \ d(1,2) \ d(1,3) \ vdots \ d(20,20) end{bmatrix}
    = begin{bmatrix} text{PSF}_{alpha,beta}(1-r_0,1-c_0) & 1 \ text{PSF}_{alpha,beta}(1-r_0,2-c_0) & 1 \ text{PSF}_{alpha,beta}(1-r_0,3-c_0) & 1 \ vdots & vdots \ text{PSF}_{alpha,beta}(20-r_0,20-c_0) & 1 end{bmatrix} times begin{bmatrix}a \ b end{bmatrix}
    + begin{bmatrix}epsilon(1,1) \ epsilon(1,2) \ epsilon(1,3) \ vdots \ epsilon(20,20) end{bmatrix}$$



    So we can write for the 1D data vector "d" :



    $$d=H(nu),theta + epsilon$$ with $H$ the matrix defined above.



    I know that we have the relation for posterior function ($d$=data and $p$ = vector of parameters :



    $$f(p|d) = dfrac{f(p),f(d|p)}{int_{p}f(d|p),text{d}p}$$



    We can also write :



    $$f(p|d) propto text{Likelihood(d|p)},f(p)quad(1)$$ with $f(p)$ the prior function (that I can take as uniform distribution).



    Now, how can I estimate the parameters $p$ from this relation $(1)$ with MCMC methods, especially in my case with Metropolis algorithm ?



    With likelihood method, I used previously the following cost function to estimate these 6 parameters :



    function cost = Crit_J(p,D)
    

    
% Compute the model corresponding to parameters p                                        
    
[R,C] = size(D);
    
[Cols,Rows] = meshgrid(1:C,1:R);                                                         
    
% Model
    
Model = (1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6));                             
    
model = Model(:);                                                                        
    
d = D(:);
    
% Introduce H matrix 
    
H = [ model, ones(length(model),1)];                                                     
    
% Compute the cost function : taking absolute value                                      
    
cost = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));                                      
    

    
end


    And after, I perform these estimations with "Matlab fminsearch" function to find a local minimum.



    Here, at first sight, I thought that I have to compute the distribution for each $p$ parameters generated randomly, but it seems to be more subtl.



    Question 1) How can we proove that, starting from arbitrary values for parameters $p$, the Metropolis algorithm will converge to right estimations ?



    Question 2) How to implement it for this concrete example ?



    Any help is welcome, regards






    maximum-likelihood monte-carlo log-likelihood







    share|cite|improve this question
























      up vote
      0
      down vote

      favorite









      up vote
      0
      down vote

      favorite











      I have 6 parameters to estimate : $p=(theta=[a,b]$, $nu=[r_0,c_0,alpha,beta])$ with Bayesian and MCMC methods :



      $$text{PSF}(r,c) = bigg(1 + dfrac{r^2 + c^2}{alpha^2}bigg)^{-beta}$$



      and the modeling :



      $$d(r,c) = a cdot text {PSF}_{alpha,beta} (r-r_0,c-c_0) + b + epsilon (r, c)$$



      with $epsilon$ being a white gaussian noise.



      This is the matrix form :



      $$begin{bmatrix}d(1,1) \ d(1,2) \ d(1,3) \ vdots \ d(20,20) end{bmatrix}
      = begin{bmatrix} text{PSF}_{alpha,beta}(1-r_0,1-c_0) & 1 \ text{PSF}_{alpha,beta}(1-r_0,2-c_0) & 1 \ text{PSF}_{alpha,beta}(1-r_0,3-c_0) & 1 \ vdots & vdots \ text{PSF}_{alpha,beta}(20-r_0,20-c_0) & 1 end{bmatrix} times begin{bmatrix}a \ b end{bmatrix}
      + begin{bmatrix}epsilon(1,1) \ epsilon(1,2) \ epsilon(1,3) \ vdots \ epsilon(20,20) end{bmatrix}$$



      So we can write for the 1D data vector "d" :



      $$d=H(nu),theta + epsilon$$ with $H$ the matrix defined above.



      I know that we have the relation for posterior function ($d$=data and $p$ = vector of parameters :



      $$f(p|d) = dfrac{f(p),f(d|p)}{int_{p}f(d|p),text{d}p}$$



      We can also write :



      $$f(p|d) propto text{Likelihood(d|p)},f(p)quad(1)$$ with $f(p)$ the prior function (that I can take as uniform distribution).



      Now, how can I estimate the parameters $p$ from this relation $(1)$ with MCMC methods, especially in my case with Metropolis algorithm ?



      With likelihood method, I used previously the following cost function to estimate these 6 parameters :



      function cost = Crit_J(p,D)
      

      
% Compute the model corresponding to parameters p                                        
      
[R,C] = size(D);
      
[Cols,Rows] = meshgrid(1:C,1:R);                                                         
      
% Model
      
Model = (1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6));                             
      
model = Model(:);                                                                        
      
d = D(:);
      
% Introduce H matrix 
      
H = [ model, ones(length(model),1)];                                                     
      
% Compute the cost function : taking absolute value                                      
      
cost = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));                                      
      

      
end


      And after, I perform these estimations with "Matlab fminsearch" function to find a local minimum.



      Here, at first sight, I thought that I have to compute the distribution for each $p$ parameters generated randomly, but it seems to be more subtl.



      Question 1) How can we proove that, starting from arbitrary values for parameters $p$, the Metropolis algorithm will converge to right estimations ?



      Question 2) How to implement it for this concrete example ?



      Any help is welcome, regards






      maximum-likelihood monte-carlo log-likelihood







      share|cite|improve this question













      I have 6 parameters to estimate : $p=(theta=[a,b]$, $nu=[r_0,c_0,alpha,beta])$ with Bayesian and MCMC methods :



      $$text{PSF}(r,c) = bigg(1 + dfrac{r^2 + c^2}{alpha^2}bigg)^{-beta}$$



      and the modeling :



      $$d(r,c) = a cdot text {PSF}_{alpha,beta} (r-r_0,c-c_0) + b + epsilon (r, c)$$



      with $epsilon$ being a white gaussian noise.



      This is the matrix form :



      $$begin{bmatrix}d(1,1) \ d(1,2) \ d(1,3) \ vdots \ d(20,20) end{bmatrix}
      = begin{bmatrix} text{PSF}_{alpha,beta}(1-r_0,1-c_0) & 1 \ text{PSF}_{alpha,beta}(1-r_0,2-c_0) & 1 \ text{PSF}_{alpha,beta}(1-r_0,3-c_0) & 1 \ vdots & vdots \ text{PSF}_{alpha,beta}(20-r_0,20-c_0) & 1 end{bmatrix} times begin{bmatrix}a \ b end{bmatrix}
      + begin{bmatrix}epsilon(1,1) \ epsilon(1,2) \ epsilon(1,3) \ vdots \ epsilon(20,20) end{bmatrix}$$



      So we can write for the 1D data vector "d" :



      $$d=H(nu),theta + epsilon$$ with $H$ the matrix defined above.



      I know that we have the relation for posterior function ($d$=data and $p$ = vector of parameters :



      $$f(p|d) = dfrac{f(p),f(d|p)}{int_{p}f(d|p),text{d}p}$$



      We can also write :



      $$f(p|d) propto text{Likelihood(d|p)},f(p)quad(1)$$ with $f(p)$ the prior function (that I can take as uniform distribution).



      Now, how can I estimate the parameters $p$ from this relation $(1)$ with MCMC methods, especially in my case with Metropolis algorithm ?



      With likelihood method, I used previously the following cost function to estimate these 6 parameters :



      function cost = Crit_J(p,D)
      

      
% Compute the model corresponding to parameters p                                        
      
[R,C] = size(D);
      
[Cols,Rows] = meshgrid(1:C,1:R);                                                         
      
% Model
      
Model = (1+((Rows-p(3)).^2+(Cols-p(4)).^2)/p(5)^2).^(-p(6));                             
      
model = Model(:);                                                                        
      
d = D(:);
      
% Introduce H matrix 
      
H = [ model, ones(length(model),1)];                                                     
      
% Compute the cost function : taking absolute value                                      
      
cost = abs((d-H*[p(1),p(2)]')'*(d-H*[p(1),p(2)]'));                                      
      

      
end


      And after, I perform these estimations with "Matlab fminsearch" function to find a local minimum.



      Here, at first sight, I thought that I have to compute the distribution for each $p$ parameters generated randomly, but it seems to be more subtl.



      Question 1) How can we proove that, starting from arbitrary values for parameters $p$, the Metropolis algorithm will converge to right estimations ?



      Question 2) How to implement it for this concrete example ?



      Any help is welcome, regards






      maximum-likelihood monte-carlo log-likelihood




      maximum-likelihood monte-carlo log-likelihood






      share|cite|improve this question













      share|cite|improve this question











      share|cite|improve this question




      share|cite|improve this question










      asked 20 hours ago









      youpilat13

      2811




      2811



























          active

          oldest

          votes











          Your Answer





          StackExchange.ifUsing("editor", function () {
          return StackExchange.using("mathjaxEditing", function () {
          StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
          StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
          });
          });
          }, "mathjax-editing");

          StackExchange.ready(function() {
          var channelOptions = {
          tags: "".split(" "),
          id: "69"
          };
          initTagRenderer("".split(" "), "".split(" "), channelOptions);

          StackExchange.using("externalEditor", function() {
          // Have to fire editor after snippets, if snippets enabled
          if (StackExchange.settings.snippets.snippetsEnabled) {
          StackExchange.using("snippets", function() {
          createEditor();
          });
          }
          else {
          createEditor();
          }
          });

          function createEditor() {
          StackExchange.prepareEditor({
          heartbeatType: 'answer',
          convertImagesToLinks: true,
          noModals: true,
          showLowRepImageUploadWarning: true,
          reputationToPostImages: 10,
          bindNavPrevention: true,
          postfix: "",
          imageUploader: {
          brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
          contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
          allowUrls: true
          },
          noCode: true, onDemand: true,
          discardSelector: ".discard-answer"
          ,immediatelyShowMarkdownHelp:true
          });


          }
          });














           

          draft saved


          draft discarded


















          StackExchange.ready(
          function () {
          StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f3004781%2fmonte-carlo-markov-chain-metropolis-hastings-estimation-of-parameters%23new-answer', 'question_page');
          }
          );

          Post as a guest















          Required, but never shown






























          active

          oldest

          votes













          active

          oldest

          votes









          active

          oldest

          votes






          active

          oldest

          votes
















           

          draft saved


          draft discarded



















































           


          draft saved


          draft discarded














          StackExchange.ready(
          function () {
          StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f3004781%2fmonte-carlo-markov-chain-metropolis-hastings-estimation-of-parameters%23new-answer', 'question_page');
          }
          );

          Post as a guest















          Required, but never shown





















































          Required, but never shown














          Required, but never shown












          Required, but never shown







          Required, but never shown

































          Required, but never shown














          Required, but never shown












          Required, but never shown







          Required, but never shown







          Popular posts from this blog

          'app-layout' is not a known element: how to share Component with different Modules

          Image
          0 My application was working fine until I decided to go with lazy loading: So, my shared component looks like this: import { Component, Renderer2 } from '@angular/core'; export interface FormModel { captcha?: string; } @Component({ selector: 'app-layout', templateUrl: './app-layout.component.html', styleUrls: ['./app-layout.css'] }) export class FullLayoutComponent { } app-layout.component.ts I am using this component in my Another component say school-home.component.html like this: <app-layout> </app-layout> and it was working fine until I created a Module school-home.module.ts like this: import { NgModule } from '@angular/core'; import { CommonModule } from '@angular/common'; import { SchoolhomeRoutingModule ...
          Read more

          android studio warns about leanback feature tag usage required on manifest while using Unity exported app?

          Image
          .everyoneloves__top-leaderboard:empty,.everyoneloves__mid-leaderboard:empty,.everyoneloves__bot-mid-leaderboard:empty{ height:90px;width:728px;box-sizing:border-box; } 0 This is a simple project showing cubes on Unity, exported to android project and imported on Android studio. On maifest file it showed following warning: Expecting <uses-feature android:name="android.software.leanback" android:required="false" /> tag I have no information about it. And also on default emptyview android project does not use/mention this feature. I know it is required by TV's and such but why its not added by default? Why is this? And Should I concern about it? android unity3d ...
          Read more

          WPF add header to Image with URL pettitions [duplicate]

          Image
          up vote 0 down vote favorite This question already has an answer here: Image source URI requiring authorization 3 answers Add header to all* HTTP requests on a machine 1 answer first of all sorry for my english. Recently I added some security in the backend of my project, and now I need to send a validation code(apiKey) inside the headers to retrieve the images. I have a problem with ALL images in my project that I had attached in my UI. For example, before I had: <Image Height="90" Source="{Binding Video.Thumbnail}"/> Now I need to change my code to create a CUSTOM petition add the apiKey in the header:...
          Read more