# Difference between revisions of "stat340s13"

m (→Multiplicative Congruential Algorithm) |
(→Class 2 - Thursday, May 9) |
||

Line 62: | Line 62: | ||

|} | |} | ||

− | ==Class 2 - Thursday, May 9== | + | ==Sampling (Generating random numbers), Class 2 - Thursday, May 9== |

− | === | + | === Introduction === |

Some people believe that activities such as rolling a dice and flipping a coin are not truly random but are '''deterministic''' – that is, the result can be reliably calculated using things such as physics and math. | Some people believe that activities such as rolling a dice and flipping a coin are not truly random but are '''deterministic''' – that is, the result can be reliably calculated using things such as physics and math. | ||

Line 73: | Line 73: | ||

=== Multiplicative Congruential Algorithm === | === Multiplicative Congruential Algorithm === | ||

− | This is an algorithm used to generate uniform, pseudo random numbers. | + | This is an algorithm used to generate uniform, pseudo random numbers. It is also referred to as the Linear or Mixed Congruential Methods. The Multiplicative Congruential Method may also refer to the special case where c=0.<br /> |

− | It is also referred to as the Linear or Mixed Congruential Methods. | ||

− | The Multiplicative Congruential Method may also refer to the special case where c=0.<br /> | ||

Line 85: | Line 83: | ||

x = x mod m – we take the remainder and go remainder mod m, but we will get x every time we do this<br /> | x = x mod m – we take the remainder and go remainder mod m, but we will get x every time we do this<br /> | ||

We can modify this:<br /> | We can modify this:<br /> | ||

− | Z = ax+b mod m<br /> | + | Z = (ax+b) mod m<br /><br /> |

− | Example: a=2, b=1, m=3; if x=10<br /> | + | |

− | Step 1: 0 = 2(10)+1 mod 3<br /> | + | '''Example''': a=2, b=1, m=3; if x=10<br /> |

− | Step 2: 1 = 2(0) + 1 mod 3<br /> | + | Step 1: 0 = (2(10)+1) mod 3<br /> |

− | Step 3: 0 = 2(1) + 1 mod 3<br /> | + | Step 2: 1 = (2(0)+1) mod 3<br /> |

+ | Step 3: 0 = (2(1)+1) mod 3<br /> | ||

You will get a sequence of numbers<br /> | You will get a sequence of numbers<br /> | ||

− | + | '''MatLab for Multiplicative Congruential Algorithm:'''<br /> | |

Before your start:<br /> | Before your start:<br /> | ||

Clear all<br /> | Clear all<br /> | ||

Close all<br /> | Close all<br /> | ||

− | |||

A=17<br /> | A=17<br /> | ||

Line 105: | Line 103: | ||

Mod(a*x+b,m)<br /> | Mod(a*x+b,m)<br /> | ||

Ans = 26<br /> | Ans = 26<br /> | ||

− | X=mod(a*x+b,m)<br /> | + | X=mod(a*x+b,m)<br /><br /> |

− | |||

+ | ''(Note: Keep repeating this command over and over again and you will seem to get random numbers – this is how the command rand works in a computer.)<br /><br /> | ||

X=mod(ax+b,m)<br /> | X=mod(ax+b,m)<br /> | ||

Line 113: | Line 111: | ||

X(ii)=mod(a*x(ii-1)+b,m);<br /> | X(ii)=mod(a*x(ii-1)+b,m);<br /> | ||

end<br /> | end<br /> | ||

− | Note: This generates 1000 numbers. <br /> | + | ''(Note: This generates 1000 numbers.)'' <br /> |

Size(x)<br /> | Size(x)<br /> | ||

Hist(x)<br /> | Hist(x)<br /> | ||

Line 127: | Line 125: | ||

− | This algorithm involves three integer parameters a, b, and m and an initial value, x_0 called seed. A sequence of numbers is defined by x(k+1) = a*x(k) + b mod m | + | This algorithm involves three integer parameters a, b, and m and an initial value, x_0 called seed. A sequence of numbers is defined by x(k+1) = a*x(k) + b mod m. Mod m means take the remainder after division by m. |

− | Mod m means take the remainder after division by m. | ||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

+ | '''Example''': a=13, b=0, m=31<br /> | ||

+ | The first 30 numbers in the sequence are a permutation of integers for 1 to 30 and then the sequence repeats itself. Values are between 0 and m-1. If the values are normalized by dividing by m-1, then result is numbers uniformly distributed in the interval [0,1]. There are only a finite number of values (30 in this case). <br /> | ||

For many years “rand” function in Matlab used this algorithm with these parameters | For many years “rand” function in Matlab used this algorithm with these parameters | ||

Line 145: | Line 136: | ||

=== Inverse Transform Method === | === Inverse Transform Method === | ||

− | U ~ U[0,1] | + | U ~ U[0,1]<br /> |

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | |||

− | x=- | + | '''Theorem''': <br /> |

+ | Take U ~ U[0,1] and let x=F^(-1)(U)<br /> | ||

+ | Then x has distribution function F(.)<br /> | ||

+ | Where F(x) = P(X<x) cdf; F^(-1)(U) denotes the function inverse of F(.)<br /> | ||

− | + | ''Ie. F(x)=U -> x=F^(-1)(U)<br /> | |

− | F^(- | + | '''Example''': f(x) = λe^(-λx)<br /> |

+ | F(x)=∫0tox〖f(x)dx〗<br /> | ||

+ | = ∫0tox〖λe^(-λx) dx〗<br /> | ||

+ | =λ/(-λ) e^(-λx) [0 to x]<br /> | ||

+ | =(-e^(-λx)+e^0 )<br /> | ||

+ | =1-e^(-λx)<br /> | ||

− | |||

− | + | y=1-e^(-λx);<br /> | |

+ | 1-y=e^(-λx);<br /> | ||

+ | x=-ln(1-y)/λ;<br /> | ||

+ | y=-ln(1-x)/λ;<br /> | ||

+ | F^(-1) (x)=-ln(1-x)/λ;<br /> | ||

+ | Step 1: Draw U ~U[0,1];<br /> | ||

+ | Step 2: x=-ln(1-U)/ λ;<br /> | ||

− | MatLab: | + | '''MatLab for Inverse Transform Method''':<br /> |

− | + | ''Clear all''<br /> | |

+ | ''Close all''<br /> | ||

− | rand(1,1000) | + | rand<br /> |

+ | rand(1,1000)<br /> | ||

+ | u=rand(1,1000)<br /> | ||

− | u | + | hist(u)<br /> |

− | + | F(x)=P(X<=x)<br /> | |

+ | =P(F^(-1)(U)<=x)<br /> | ||

+ | =P(F(F^(-1)(U))<=F(x))<br /> | ||

+ | =P(U<=F(x))<br /> | ||

− | + | <!-- Did class end before this was finished? --> | |

− | |||

− | |||

− |

## Revision as of 12:47, 9 May 2013

## Contents

**Computer Simulation of Complex Systems (Stat 340 - Spring 2013)**

## Introduction, Class 1 - Tuesday, May 7

### Four Fundamental Problems

1. Classification: Given an input object X, we have a function which will take in this input X and identity which 'class (Y)' it belongs to (Discrete Case)

2. Regression: Same as classification but in the continuous case

3. Clustering: Use common features of objects in same class or group to form clusters.

4. Dimensionality Reduction

### Applications

Most useful when structure of the task is not well understood but can be characterized by a dataset with strong statistical regularity

Examples:

- Computer Vision, Computer Graphics, Finance (fraud detection), Machine Learning
- Search and recommendation (eg. Google)
- Automatic speech recognition, speaker verification
- Text parsing
- Face identification
- Tracking objects in video
- Financial prediction, fraud detection

### Course Information

**General Information**

- No required textbook, recommended: "Simulation" by Sheldon M. Ross
- Computing parts of the course will be done in Matlab, but prior knowledge of Matlab is not essential (will have a tutorial on it)
- Learn for announcements, assignments, and emails.
- Other course material on: http://wikicoursenote.com/wiki/
- Log on to both Learn and wikicoursenote frequently.

**Primary contributor**: Put a summary of the lecture up within 48 hours.

**General contributor**: Elaborate on concepts, add example, add code, add pictures, reference, corrections etc… withing 2 weeks

~~Must do both~~ *All primary contributions are now considered general contributions you must contribute to 50% of lectures for full marks*

- A general contribution can be correctional (fixing mistakes) or technical (expanding content, adding examples, etc) but at least half of your contributions should be technical for full marks

Do not submit copyrighted work without permission, cite original sources.
Each time you make a contribution, check mark the table. Marks are calculated on honour system, although there will be random verifications. If you are caught claiming to contribute but didn't, you will *lose* marks.

### Tentative Marking Scheme

Item | Value |
---|---|

Assignments (~6) | 30% |

WikiCourseNote | 10% |

Midterm | 20% |

Final | 40% |

## Sampling (Generating random numbers), Class 2 - Thursday, May 9

### Introduction

Some people believe that activities such as rolling a dice and flipping a coin are not truly random but are **deterministic** – that is, the result can be reliably calculated using things such as physics and math.

A computer cannot generate truly random numbers since computers can only run algorithms, which are deterministic in nature. They can, however, generate **Pseudo Random Numbers**; numbers that seem random but are actually deterministic.

### Multiplicative Congruential Algorithm

This is an algorithm used to generate uniform, pseudo random numbers. It is also referred to as the Linear or Mixed Congruential Methods. The Multiplicative Congruential Method may also refer to the special case where c=0.

Take a number x and divide it by m

The remainder is a number between 0 and m-1

We use the operator “mod”

X mod m

1 = 10 mod 3

x = x mod m – we take the remainder and go remainder mod m, but we will get x every time we do this

We can modify this:

Z = (ax+b) mod m

**Example**: a=2, b=1, m=3; if x=10

Step 1: 0 = (2(10)+1) mod 3

Step 2: 1 = (2(0)+1) mod 3

Step 3: 0 = (2(1)+1) mod 3

You will get a sequence of numbers

**MatLab for Multiplicative Congruential Algorithm:**

Before your start:

Clear all

Close all

A=17

B=3

M=31

X=5

Mod(a*x+b,m)

Ans = 26

X=mod(a*x+b,m)

*(Note: Keep repeating this command over and over again and you will seem to get random numbers – this is how the command rand works in a computer.)*

X=mod(ax+b,m)

For ii=2:1000

X(ii)=mod(a*x(ii-1)+b,m);

end

*(Note: This generates 1000 numbers.)*

Size(x)

Hist(x)

X(1)=5

For ii=2:5000

X(ii)=mod(a*x(ii-1)+b,m);

End

Hist(x)

*(Note: The semicolon after the x(ii)=mod(a*x(ii-1)+b,m) ensures that Matlab will not show the entire vector that is x. It will instead calculate it internally and you will be able to work with it. Adding the semicolon to the end of this line reduces the run time significantly.) *

This algorithm involves three integer parameters a, b, and m and an initial value, x_0 called seed. A sequence of numbers is defined by x(k+1) = a*x(k) + b mod m. Mod m means take the remainder after division by m.

**Example**: a=13, b=0, m=31

The first 30 numbers in the sequence are a permutation of integers for 1 to 30 and then the sequence repeats itself. Values are between 0 and m-1. If the values are normalized by dividing by m-1, then result is numbers uniformly distributed in the interval [0,1]. There are only a finite number of values (30 in this case).

For many years “rand” function in Matlab used this algorithm with these parameters A=7^5=16807, b=0, m=2^31 -1, 2147483647 – recommended in a 1988 paper by Park and Miller (Important part is that m should be large)

### Inverse Transform Method

U ~ U[0,1]

**Theorem**:

Take U ~ U[0,1] and let x=F^(-1)(U)

Then x has distribution function F(.)

Where F(x) = P(X<x) cdf; F^(-1)(U) denotes the function inverse of F(.)

*Ie. F(x)=U -> x=F^(-1)(U)*

**Example**: f(x) = λe^(-λx)

F(x)=∫0tox〖f(x)dx〗

= ∫0tox〖λe^(-λx) dx〗

=λ/(-λ) e^(-λx) [0 to x]

=(-e^(-λx)+e^0 )

=1-e^(-λx)

y=1-e^(-λx);

1-y=e^(-λx);

x=-ln(1-y)/λ;

y=-ln(1-x)/λ;

F^(-1) (x)=-ln(1-x)/λ;

Step 1: Draw U ~U[0,1];

Step 2: x=-ln(1-U)/ λ;

**MatLab for Inverse Transform Method**:

*Clear all*

*Close all*

rand

rand(1,1000)

u=rand(1,1000)

hist(u)

F(x)=P(X<=x)

=P(F^(-1)(U)<=x)

=P(F(F^(-1)(U))<=F(x))

=P(U<=F(x))