Info

La pregunta está cerrada. Vuélvala a abrir para editarla o responderla.

Hi could someone help me figure out how to use this script for MonteCarlo please?

1 visualización (últimos 30 días)
Michael Lee
Michael Lee el 14 de Mayo de 2017
Cerrada: John D'Errico el 14 de Mayo de 2017
%H62NUM - Numerical Methods and Contextual Electrical and Electronic
%Engineering Mathematics
%Dr A J Phillips March 2016
%Script file that performs Monte Carlo simulation
clc %to clear the Command Window before beginning so everything we see relates to the latest run
clear all %to clear all previously existing values in memory so the outcomes of previous scripts or commands do not impact on the running of this script
%EDITABLE CODE - You MAY edit the number of trials and the number of bins
trials=100000; %the overall number of trials ("experiments"). More is better in terms of accuracy but not in terms of speed!
nobins=100; %how many points do we want on our final pdf (how many bars on our histogram)?
%EDITABLE CODE - You SHOULD set any parameters relevant for your chosen random variables here, BEFORE
%the trials loop. This applies to both the case where we use a formula to
%sample the pdf (via the cdf) and where we need to generate a vector
%representation of the cdf for interpolation.
%lambda= 1; %a possible parameter for an exponential (negative exponential) distribution
%sigma=; %a possible parameter (standard deviation) for a Gaussian (normal) distribution
%mu=; %a possible parameter (mean) for a Gaussian (normal) distribution
%a=; %a possible lower value for a scaled and shifted uniform distribution
%b=; %a possible upper value for a scaled and shifted uniform distribution
a=1;
b=1;
sigma=1;
mu = 0;
lambda= 1;
%If any random variable needs a cdf generating do it here
%EDITABLE CODE - You WILL in general need to change x1low and x1high (and similar
%values for further random variables)
%x1low=; %the minimum value x1 can take (either the true minimum based on there being an actual lower limit or a constraint placed on a very long thin and possibly infinite pdf lower tail)
%x1high=; %the maximum value x1 can take (either the true maximum based on there being an actual upper limit or a constraint placed on a very long thin and possibly infinite pdf upper tail)
%if too wide a range is chosen however then the computer will assign cdf values
%of precisely 1 for high x1 and precisely x1high=0 for low x1 and will be unable
%to interpolate later
%EDITABLE CODE - You MAY wish to change the number of steps used in the cdf
%integration
%steps=10000; %controls the spacing (delx1) used for the initial integration of the pdf to get the cdf. It is only done once (per randome variable) and so it is ok to use a decent number of steps.
%DO NOT EDIT delx1 and x1vec (except to delx2=(x2high-x2low)/steps; x2vec=x2low:delx2:x2high; etc.)
%delx1=(x1high-x1low)/steps;
%x1vec=x1low:delx1:x1high;
%EDITABLE CODE - You WILL need to put here any pdf that you need to
%integrate numerically to get the cdf
%pdfin1vec=(exp(-((x1vec-mu).^2)/(2*sigma^2)))/(sqrt(2*pi)*sigma);
%DO NOT EDIT cdf1vec (except to cdf2vec=cdfcalc(x2vec,pdfin2vec);etc.
%cdf1vec=cdfcalc(x1vec,pdfin1vec);
% x1low=-5;
% x1high=5;
% steps=1000;
% delx1=(x1high-x1low)/steps;
% x1vec=x1low:delx1:x1high;
% pdfin1vec=(exp(-((x1vec-mu).^2)/(2*sigma^2)))/(sqrt(2*pi)*sigma);
% cdf1vec=cdfcalc(x1vec,pdfin1vec);
%EDITABLE CODE
%given both the individual pdfs you are sampling and how you are operating
%and combining the random variables you should now be able to work out the
%range of possible outcomes
low=-6; %the minimum value our overall randomtrialval can take (either the true minimum based on there being an actual lower limit or a constraint placed on a very long thin and possibly infinite pdf lower tail)
high=6; %the maximum value our overall randomtrialval can take (either the true maximum based on there being an actual upper limit or a constraint placed on a very long thin and possibly infinite pdf upper tail)
%in setting low and high of course need to consider x1low, x1high (etc.)
%and any transformation or combination of the random variables (and be
%particularly careful to consider the range of random variables for whom we
%have gone straight to use of a formula in the loop below
%DO NOT EDIT
binwidth=(high-low)/nobins;
bin1centre=low+binwidth/2;
binvec=1:nobins;
bincentrevec=bin1centre:binwidth:bin1centre+(nobins-1)*binwidth;
count=zeros(1,nobins); %creates an array of all zeros, essentially initializing count(1)=0, count(2)=0,..., count(nobins)=0
count_out_of_range=0; %sometimes our bins cannot quite cover the whole range of a random variable (see definition of low and high above)
%EDITABLE CODE - You WILL need to uncomment sumforav, sumfortailprob and threshold at some stage.
%sumforav=0; %for expectation calculation direct from the trials loop (as opposed to from the eventual histogram/pdf)
%sumfortailprob=0; %for tail probability calculation direct from the trials loop (as opposed to from the eventual histogram/pdf)
%threshold=; %for tail probability
for i=1:trials
%EDITABLE CODE - You WILL need to edit below to include,
%combine and remove various input random variables.
%The basic format of u1=rand; followed by either a formula derived from the
%cdf or interpolation directly using the cdf must be retained for each
%random variable used.
%Random variables requiring interpolation from a cdf vector
%u1=rand; %uniform random number between 0 and 1 generated here
%x1=interp1(cdf1vec,x1vec,u1);
%etc. for more random variables
%u2=rand; %uniform random number between 0 and 1 generated here
%x2=interp1(cdf2vec,x2vec,u2);
%If you get an error about "The grid vectors are not strictly monotonic
%increasing" then look at the definition of low and high
%Random variables permitting the use of a formula
%u3=rand; %uniform random number between 0 and 1 generated here
%x3=; %rhs is function of u3
%EXAMPLE (1) x3=-log(1-u3)/lambda works for the exponential pdf p(x)=lambda*exp(-lambda*x) x>=0 which has cdf F(x)=1-exp(-lambda*x)
%EXAMPLE (2) x3=a+u3*(b-a) works for the scaled uniform between a and b (a < b) pdf p(x)=1/(b-a) a<x<b which has cdf F(x)=(x-a)/(b-a)
%etc. for more random variables
%u4=rand; %uniform random number between 0 and 1 generated here
%x4=; %rhs is function of u4
u1=rand; %generates a uniform random variable on (0,1)
x1=u1*(b-a);
u2=rand; %generates a uniform random variable on (0,1)
x2=u2*(b-a);
x3 = x1+x2;
%in general a transformation from a uniform random variable on (0,1)
%x1=interp1(cdf1vec,x1vec,u1);
randomtrialval=x1; %this may be a function of x1, x2, x3 etc. (f(x1,x2,x3...) [don't use the f notation, just type in 2*x1+x2^2 or whatever function you want]
%sumforav=sumforav+randomtrialval; %this is preparing us to find E(f(x1,x2,x3...)
% if randomtrialval>threshold %finding the upper tail probability (use < for a lower tail)
% sumfortailprob=sumfortailprob+1;
% end
%finding the appropriate bin whose count we then increment
%DO NOT EDIT this part of the loop
binno=round(((randomtrialval-bin1centre)/binwidth)+1);
if ( (binno>=1) && (binno<=nobins))
count(binno)=count(binno)+1;
else
count_out_of_range=count_out_of_range+1;
end
end
%av=sumforav/trials;
%tailprob=sumfortailprob/trials;
pdfvec=count/(trials*binwidth);
figure(1)
plot(bincentrevec,pdfvec);
xlabel('Random trial value')
ylabel('pdf')
%EDITABLE CODE: If you (optionally) want to see the original histogram (as a line plot) uncomment the code below
figure(2)
plot(binvec,count)
xlabel('Bin number')
ylabel('Bin count')
prob_out_of_range=count_out_of_range/trials %if a substantial fraction of 1 you should reconsider your bin definition
sumarea=0;
for i=1:nobins
sumarea=sumarea+binwidth*pdfvec(i); %simply a quick check that the pdf does have an area of 1.
%As it is a quick check we accept the poorer accuracy provided by the midpoint integration approach (i.e. we may not always get *precisely* one here).
end
pdfareacheck=sumarea %Should be equal to 1
function [F] = cdfcalc(x,pdf )
%DO NOT EDIT THIS FUNCTION
%integrate the pdf to get the cdf in very small strips using trapezoidal rule for simplicity
F(1)=0; %no integration is performed for the first element in the pdf.
%If this has been unwisely chosen (i.e where there is an actual significant negative tail before it) it will constitute an innaccuracy.
sum=0;
for i=2:length(x)
striparea=(x(i)-x(i-1))*(pdf(i)+pdf(i-1))/2;
sum=sum+striparea;
F(i)=sum;
end
end
Need help simulating the tasks in the image, please help!

Respuestas (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by