COMPUTL BIOL & CHEM
COMPUTL BIOL & CHEM AMATH 410
Popular in Course
Popular in Applied Mathematics
This 16 page Class Notes was uploaded by Mossie Pfannerstill V on Wednesday September 9, 2015. The Class Notes belongs to AMATH 410 at University of Washington taught by Staff in Fall. Since its upload, it has received 9 views. For similar materials see /class/192278/amath-410-university-of-washington in Applied Mathematics at University of Washington.
Reviews for COMPUTL BIOL & CHEM
Report this Material
What is Karma?
Karma is the currency of StudySoup.
Date Created: 09/09/15
AMATH 410 Review Structuring a program in MATLAB Our aim is to illustrate a basic strategy for how to write a new MATLAB program from scratch A m le with the code below is available as cellreproductionreviewcode1m from our website H de ne given parameters 2 make a dummy list or matrix of values for the variable that are going to evolve over time 3 set the initial value of that dummy list 4 run a for loop to ll in all values of that variable solution list 5 process the result7 if required 6 plot Lets say were given this problem 7 similar to the cell reproduction model discussed on the rst day of class Implement the rule nt 4pnt7117 nt 7 starting from an initial population size of 711 10007 and simulating for 100 timesteps or generations Then7 compute a list of all times t where nt gt 1300 Use p 0757 K 2000 Step 1 Where do I start In my rst lines of code7 I de ne all of the parameters p75 oofraction of cells that survive each reproduction generationend100 oopredict up to this generation in the future K2000 ooreproduction capacity Step 2 OK7 now I have to think about what it is that I am going to evolve over time I am going to prede ne a dummy77 list or matrix of numbers representing the variable that I am going to simulate That dummy list is often just lled in with zeros Then I will go through and ll it in later Here7 I want a list of values of n at timesteps t 17 100 So7 lets just de ne a list of 100 zeros to use as my dummy list There are a couple of ways to do this nlistzeros 1 generationend or nlist0 1 generationend Step 3 Alright Now I need to set my initial conditions Whats the rst value of nlist That7s nlist11000 Step 4 Now I need to loop over time At each timestep7 l implement my dynamical rule for t2zgenerationend nlistt4nlistt 1p1 nlistt 1K end So7 I went through and lled in the dummy zero values of nlist with their real values Step 5 Now7 I need to make a new list 7 of generations t when nlist tgt30 Let7s call this new list listoftimesoverlimit Note that I DONT KNOW IN ADVANCE HOW LONG THIS LlST lS SUPPOSED TO BE In this way7 the situation is similar to when we are building up a list of dwell times for the markov chain We will therefore not make a dummy list in advance7 but proceed as follows all three ways give same answer Way 1 i keep growing the list Via a loop Here7 we have to de ne an integer valued number whichentrytoupdatem that tells me which entry in the list l7m supposed to be changing We start with the rst entry whichentrytoupdate1 OK7 let7s use this for t1generationend if nlisttgt1300 listoftimesoverlimitwhichentrytoupdatet oonext I say update the NEXT entry down the list next time whiChentrytoupdatewhiChentrytoupdate1 end end Way 2 i keep growing the list Via a loop version 2 Here7 we start with an empty list7 and keep appending to it listoftimesoverlimit for t1generationend if nlisttgt1300 listoftimesoverlimitlistoftimesoverlimit t end end Way 3 i nd command Here I use the find command in matlab listoftimesoverlimitfindnlistgt1300 That7s all we have to do Here7s a review on using the nd command find mylist3 ooreturns a list of indices for entries that are equal to 3 find mylistlt10 ooreturns a list of indices for entries that are less than 10 find mylist 5 ooreturns a list of indices for entries that are not equal to 5 The point is you type find CONDITION and get out all the entries where CONDITION is true For example7 find 1 2 2 2 returns 2 3 Step 6 OK7 let7s plot 71 vs t I7ll make a list of times It rst tlist1zgenerationzend figure plot tlistnlist 7 7 7MarkerSize7 24 ylabel 7nt 7 7FontSize7 20 Xlabel 7t7 7FontSize7 20 set gca 7FontSize 7 20 AMATH 410 MATLAB TECHNIQUES This is a PARTIAL list of techniques that you ll want to be familiar with If you re not please experiment in MATLAB use the help commands and ASK the course instructors in of ce hours or class lab sessions We re here to help and all of us have questions about these things In addition please make sure that you re familiar with the MATLAB functions and techniques that are listed in the sections of the Lab Manual that we have covered so far in class and the commands used in class please carefully look at all codes from class which are posted online with class notes Making m les using the text editor To create a new m le click File gt New gt m file or click the bottom right underneath File that looks like a blank sheet of paper Say you make a code called myprogramm 7 you save it in a directory make sure that this is your working directory see below and type myprogram at the command line to run this code A good habit is to write clear all close all on top of your code so that every time before it runs MATLAB clears up the workspace and existing plot windows Changing the working directory in MATLAB At the top of the matlab command window it ll say what directory you are currently working in Of course this needs to be the same one where you are saving your m program les and any data les to To change the working directory either click on the bar where this is displayed or use the cd command at the command line Use help functionname where functionname is anything you don t know how to use yet say help hist whenever you want more info doc functionname often gives even more info with examples 7 cutting and pasting these to the command line is a great way to get started If you have no idea about the function name say you are looking for a function that does something you can search in the MATLAB help document using key words whos command 7 type this at command line to see list of all variables and their sizes Or to see values of a variable either go to workspace window as described in Lab Manual or simply type its name at the command line eg type A then hit the enter key to see value of A Editing m program les just type edit and an editor will pop up to edit the code myfilename m type edit myfilename Loading data using load Commands size and length Command linspace Using the colon to extract parts of vectors or matrices Using the semi colon to hide the print outs you are not interested in viewing Using for loops Using the disp and numZStr commands to output data to screen in organized way 7 see p15 of Lab Manual Using and writing MATLAB functions Use of anonymous functions and function handles Use of the hold on command for plotting multiple curves or points etc in a single plot 7 Syntax figure plotxlist1 ylistl hold on plotxlist2 ylist2 7 Please compare this with what happens when you leave off the hold on 7 When you no longer needs the 77hold on77 property type 77hold off7 Get familiar with the rich options MATLAB has to offer in plotting by checking out the help document of 77plot Simply type 7 help plot77 in the command window and see what MATLAB likes to tell you This is also helpful when you forget which letter stands for which color or what kind of lines in plotting Use of the percent sign 7 to start a comment notation that is ignored by MATLAB Matrix multiplication via and quot Use of the vectorized operations quot etc to do ELEMENT WISE operations on entries of vectors or matrices Please note that for example the command is therefore very much different from the matrix multiplication operation of linear algebra Using format control to specify numbers of signif digits etc And in week 3 and later weeks as this material is covered Using eigenvalue eigenvector commands via eig Use of find command to identify extract or modify target entries in a matrix or vector eg sec 8 of lab manual Using max and abs 7 eg to nd dominant eigenvalues eigenvectors sec 8 of lab manual Linear solve of matrix equations Wc n where w is matrix n is given vector and 6 unknown vector that are solving for Syntax cAn Using rand and randn random number generators Using the mean and var commands Plotting and properly scaling a histogram Simulating a Markov chain 7 using logical operators and repeated calls to rand AMATH 410 Review Implementing basic probability concepts in MATLAB First7 we discuss generating samples or realizations trials of a random variable x This means using a random number generator to produce a number that occurs with the frequencies or probabilities that are pre speci ed for x These probabilities can be speci ed in one of two ways 0 First7 if x is a discrete valued random variable7 with N possible values making up the state space7 or sample space 5152775N then there is a list of N values Px 57 19 that specify the probabilities Examples of this from class are the binary and binomial random variables Second7 if z is a continuous valued random variable7 there is a whole range of possible values that it could take7 not just a discrete list well just talk about one of this type of random variable here7 the uniform distributed random variable What we mean by frequencies of occurrence77 above is7 if I generate the random variable over and over again M times for a huge number of trials M7 on what fraction of those trials will z take a certain value Let us say we count up a total of 771 times that the state 57 occurs The fact that px 57 pj means that mjM 10 That7s it 19 is the frequency or fraction of trials where z is found in a certain state ie7 where z takes the jth value7 j 8739 Counting up the total number of times that anything at all happened in our trials above7 Z 771739 M 739 This means7 by our de nition7 Z 1 739 That7s the frequency or fraction of trials when anything at all any possible state occurred That is equal to 1 That basic fact is known as normalization7 but its just logic 7 whenever we talk about the probability of anything happening7 we7d better be able to sum up the probabilities of all of the different possibilities and get 1 If this does NOT work7 then I need to x the situation One place where this comes up is in Markov chains 1 nd the eigenvector W with domi nant eigenvalue 1 Lets say there are three possible states in my markov chain7 so that Now7 the different 10 here are supposed to give us the probabilities pj of being in states 57 at equilibrium They7d better add to 1 ln general7 they will NOT add to 1 when they are produces as entries of an eigenvector7 coming straight of out of MATLAB What do I do MAKE them add to one by dividing each by the sum of the rest The corresponding probabilities are then 101 wiwi 102 ws 102 wzwi 102 ws 103 wswi 102 ws This certainly gives us p1 192 133 1 so these now make sense as probabilities The rand command The uniform distributed random variable with range 01 is the basis for most of our codes involving probability models in the course We START with it and generate samples of binary or binomial random variables Or we START with it and determine into which state a markov chain should transition at the next timestep k The uniform distributed random variable with range 01 takes values between 0 and 1 7 any value there with equal probability It is a continuous valued random variable Unlike for a discrete valued random variable its not very useful to assign a probability to a speci c value of the uniform distrubuted random variable since that requires z being equal to this value with perfect precision Eg the fraction of trials when say z 341741902437192347 is extremely small 7 actually zero if we de ne what z is supposed to be with perfect precision endless digits So instead we specify the probability that z lies anywhere in subrange ab lf0 a b 1thenwe have Px 6 ab b7a 1 To generate a single realization of such a uniformly distributed random variable with range 01 we type Xrand If we want a list of M samples we type Xrand1M If we want a M gtlt M matrix of samples we type XrandMM The rand command and binaryvalued random variables Say we want to use these samples to make a list of samples of a binary valued random variable s with Px 1 H and Px 0 1 7 H 1 start with a uniform random variable 7 this time 1711 call this 7 Then I will use equation This gives us the fact Pr E 0H So if 1 set x 1 every time that r E 0 H and z 0 otherwise x will be the desired binary valued random variable Speci cally rrand if rltH X1 else XO end does the trick Now say 1 want to make not just one sample of a binary valued random variable but many say M of them all in a list xlist Here7s the code for that see the HW solution code that solves Ex 33 rlistrand1M Xlistzeros1M oAwe just declared Xlist as a I39dummyquot or quotblankquot list oAwe re about to fill it in with real values for j1 M if rlistjltH Xlistj1 else Xlistj0 end There is also faster but absolutely equivalent way of doing this also explained in the comment annotated HW solution codes rlistrand1M Xlistzeros1M XlistfindrlistltH1 The rand command and simulating Markov chains Here we want to use rand to do something different but very closely related This is also explained in detail in the comments in the HW solution code EandGexercise36 m The situation we7ll discuss is that our Markov chain is currently in state 3 at timestep k That is 3 1 want to advance to the next timestep That means setting xk 1 1 am going to do this based on the transition probabilities that belong to the jth column of the matrix A Let us say that there are 3 states so that A is a 3 gtlt 3 matrix Then the third column is Here 113 is the probability that 1 will transition from state 3 to state 1 123 is the probability that I will transition from state 3 to state 2 133 is the probability that 1 will stay in state 3 OK so which of there three possibilities happens on a given trial Clearly I need to generate a random variable xk1 7 the state on the next timestep 7 such that zk 1 1 with probability 113 zk 1 2 with probability 123 zk 1 3 with probability 133 1 will use rand to do this Just like for the binary valued random variable above we rst set rrand Next 1 use the facts Pr E 0a13 113 Pr E a13a13 123 123 Pr E a23a23 133 133 exactly as above So splitting up the interval 01 in that way gives me the probabilities that l was after Well l7rn now going to use this to de ne the value ofthe state on the next tirnestep zk1 Here7s the setup exactly as in EandGexerCise36m l have a list called states that is we are gradually lling in with the values Say l7ve already sirnulated up to tirnestep k Then states k and we are going to set states k1xk 1 Lets do it Here we are ASSUMING that statesk3 like the above The code EandG exercise 3 6m handles the more general case where we dont prespecify the value of states k that7s a bit arti cial but simplest rrand if r lt a13 oofor transition FROM state 3 to state 1 states k11 elseif r lt a13a23 oofor transition FROM state 3 to state 2 states k12 else statesk13 oofor transition FROM state 3 to state 3 end That7s it Let me repeat what we7d do in general without the arti cial prespecifying just discussed if r lt A1statesk Zfor transition FROM statesk to state 1 statesk11 elseif r lt A1stateskA2statesk Zfor transition FROM statesk to state 2 statesk12 else statesk13 Zfor transition FROM statesk to state 3 end