modelos probabilisticos

INGENIERIA de TELECOMUNICACION

ESTADISTICA

2011-2012

PRACTICA 3. MODELOS DE PROBABILIDAD

OBJETIVOS: Introduccion a modelos de probabilidad discretos y contnuos mas comunes.Caracterizacion, representacion grafica. Resolucion mediante simulacion con MATLAB/Octave deejercicios propuestos.

1. Introduccion

En general, para generar variables aleatorias contnuas utilizaremos el metodo de la trans-formacion inversa de la funcion de distribucion (siempre que exista inversa). Para el caso de v.a.sdiscretas utilizaremos una condicion booleana. No obstante, MATLAB/Octave, dispone de funcio-nes propias para generar variables aleatorias de modelos de probabilidad conocidos.

La siguienta tabla resume algunas de las funciones mas importantes para la generacion de nume-ros aleatorios de modelos de probabilidad para variables contnuas y discretas en MATLAB/Octave:

Funcion Descripcion Sintaxis

normrnd no aleatorios N(, ) normrnd(MU,SIGMA,m,n)randn no aleatorios N(0, 1) randn(m,n)exprnd no aleatorios Exp() exprnd(1/lambda,m,n)binornd no aleatorios Bin(n, p) binornd(N,P,m,n)poissrnd no aleatorios Poiss() poissrnd(lambda,m,n)

donde m y n, son respectivamente en no de filas y columnas a generar.

2. Caso contnuo

Distribucion normal

Recordemos que la distribucion Normal tiene como funcion de densidad:

X N (, ) : f(x) = 1

2piexp

122

(x)2 , donde x R, > 0 y R.

1. Crea una funcion en MATLAB/Octave que proporcione los valores de la funcion de densidadde una variable aleatoria con distribucion normal N (, ).

Ingeniera de Telecomunicacion - Estadstica (2011-2012), PRACTICA 3. MODELOS DE PROBABILIDAD 1

%% creamos la funcion fd_normal.m

function y = fd_normal(x, mu, sigma)

for i=1:length(x)

y(i) = exp(-0.5*((x(i)-mu)/sigma)^2) / (sigma*sqrt(2*pi));

end

%%

2. Utilizando la funcion creada en el apartado anterior, representa graficamente la funcion dedensidad de la v.a. normal para distintos valores de los parametros de la distribucion:

a) Manteniendo y variando Considera 3 distribuciones Normales con desvia-cion tpica constante ( = 1) y con diferentes medias ( =-1,0,1). Como afectan losparametros a la forma de la distribucion? y a su posicion en los ejes?

b) Mantiendo y variando Considera 3 distribuciones Normales con media cons-tante ( = 0) y variando la desviacion tpica ( =0.3,0.5,1.2). Analizar como afecta a laforma de la distribucion y como afecta a su posicion en los ejes.

% a)

>> x = -5:0.01:5;

>> y1 = fd_normal(x,-1,1); % y1 es Normal con media -1 y desv. 1

>> y2 = fd_normal(x, 0,1); % y2 es Normal con media 0 y desv. 1

>> y3 = fd_normal(x, 1,1); % y2 es Normal con media +1 y desv. 1

% b)

>> y4 = fd_normal(x, 0, 0.3); % y4 es Normal con media 0 y desv. 0.3



Una vez generadas las 6 variables aleatorias normales, vamos a compararlas graficamenteusando la funciones subplot1; plot y hold on, hold off.

1mediante la funcion subplot podemos crear varios graficos a la vez. En el ejemplo, subplot(1,2,i), significa quedibujamos en 1 fila y en 2 columnas el grafico i = 1 y 2


% GRAFICO CASO a)

>> subplot(1,2,1) % 1 fila, 2 columnas, grafico no 1

>> hold on % esta opcion permite superponer graficos

>> plot(x,y1, b) % grafico x/y1 en color azul (blue)

>> plot(x,y2, g) % grafico x/y2 en color verde (green)

>> plot(x,y3, r) % grafico x/y2 en color rojo (red)

>> hold off % deshabilitamos la opcion para superponer graficos

% GRAFICO CASO b)

>> subplot(1,2,2) % 1 fila, 2 columnas, grafico no 2

>> hold on

>> plot(x,y4, b)

>> plot(x,y5, g)

>> plot(x,y6, r)

>> hold off

% Otro modo de hacerlo es:

>> subplot(1,2,1)

>> plot(x,y1,b,x,y2,g,x,y3,r)

>> subplot(1,2,2)

>> plot(x,y4,b,x,y5,g,x,y6,r)

-5 0 50

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

-3 -2 -1 0 1 2 30

0.2

0.4

0.6

0.8

1

1.2

1.4

mu = -1 , sigma = 1mu = 0 , sigma = 1mu = 1 , sigma = 1

mu = 0 ,sigma = 0.3mu = 0 , sigma = 0.5mu = 0 ,sigma = 1.2

NOTA: Otra manera de representar un grafico de una distribucion normal ( = 0 y = 2,5)para x (10, 10) sera:

>> fplot(fd_normal(x, 0, 2.5), [-10 10], g)

Con la sentencia disttool de MATLAB sobre el Command Window, podemos ver graficamente


la densidad de diferentes distribuciones, seleccionando en Function Type/PDF.

3. Genera en MATLAB/Octave numeros aleatorios de la distribucion normal.

En MATLAB/Octave es posible generar numeros aleatorios de una distribucion normal(, ), a traves de funciones propias de la librera estadstica stats. As, podemos creardatos normales (n filas y m columnas) de la siguiente manera:

>> x= normrnd(MU,SIGMA,m,n);

El comando randn, permite generar numeros aleatorios de una distribucion normal estandar(es decir, de media = 0 y = 1). A partir de randn, es posible generar numeros aleatorios N (, ). Del siguiente modo: sea Z N (0, 1), podemos obtener X N (, ), dada lasiguiente relacion:

X = Z +

>> z= randn(m,n);

>> x= z*sigma + mu;

Podemos comprobar, que los numeros generados son N (, ) representando el histograma:

>> m=1000; n=1;

>> sigma=0.75; mu=1;

>> z=randn(m,n); x=z*sigma + mu;

>> hist(x)

-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.50

50

100

150

200

250

300

4. Funcion de distribucion acumulada. En MATLAB/Octave, el comando normcdf(x,mu,sigma)devuelve la probabilidad p = P (X x) de una distribucion normal de parametros y .Representa la funcion de distribucion acumulada, para los valores de x [3, 3], siendo Xuna normal estandar.


>> x = [-3:0.01:3];

>> mu=0; sigma=1;

>> p = normcdf(x,mu,sigma); % p = normcdf(x) proporciona

% la f. distr. acum. de

% una N(0,1)

% Su inversa es x = norminv(p,mu,sigma)

>> plot(x,p)

>> grid on

-3 -2 -1 0 1 2 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

3. Caso discreto

Distribucion binomial

La funcion de probabilidad p(x) de la distribucion binomial Bin(n, p) es:

X Bin(n, p) : p(x) =(n

x

)px(1 p)nx donde x = 0, 1, ..., n y 0 p < 1.

Siendo n el numero de intentos o de ensayos y el parametro p la probabilidad de que el sucesoexito ocurra.

1. Crea una funcion en MATLAB/Octave que permita representar funciones de probabilidad deuna distribucion binomial.

%% creamos la funcion fp_binomial.m

function y=fp_binomial(x,n,p)

for i=1:length(x)

y(i)=(nchoosek(n,x(i)))*(p^x(i))*((1-p)^(n-x(i))) ;

end

%%

NOTA: la funcion nchoosek(n,x) permite calcular(nx

), sin embargo, para un tamano de

n grande, el resultado puede no ser exacto. Por ello, para un tamano n grande, emplearemos


la funcion Gamma :

(p) =

+0

exxp1dx, donde p > 0

Propiedades de la funcion Gamma ():

a) (1) = 0! = 1

b) (p) = (p 1)(p 1), p > 0 y entonces (p) = (p 1)!, p N.c)(nx

)= (n+1)(x+1)(nx+1)

d) (12) =pi.

2. Crea en MATLAB/Octave que permita representar funciones de probabilidad de una distri-bucion binomial para tamanos de n grandes.

%% creamos la funcion fp_binomialN.m

%% utilizando la propiedad 3) de la funcion Gamma

function y = fp_binomialN(x,n,p)

for i=1:length(x)

y(i)=(gamma(n+1)*(p^x(i))*((1-p)^(n-x(i))))/

(gamma(x(i)+1)*gamma(n-x(i)+1));

end

%%

3. Utilizando la funcion creada anteriormente, representa graficamente la funcion de probabilidadde una variable aleatoria binomial, para los siguientes casos:

a) Dejando constante p y variando n creamos 3 distribuciones binomiales:Bin(5,0.2), Bin(10,0.2), Bin(20,0.2)

b) Dejar n constante y variar p creamos otras 3 distribuciones binomiales:Bin(100,0.1), Bin(100,0.5), Bin(100,0.8)


% Caso a):

>> x5=0:5 % creamos una secuencia de 0 a 5

>> y5=fp_binomial(x5,5,0.2) % llamamos a la funcion fp_binomial.m

>> sum(y5) % comprobamos que la suma de las probabilidades es 1

>> x10=0:10;

>> y10=fp_binomial(x10,10,0.2);

>> x20=0:20;

>> y20=fp_binomial(x20,20,0.2);

% Caso b):

>> x100 = 0:100;

>> y1 = fp_binomialN(x100, 100, 0.1); % llamamos a la funcion

>> y2 = fp_binomialN(x100, 100, 0.5); % fp_binomialN.m

>> y3 = fp_binomialN(x100, 100, 0.8);

Podemos comparar ambos graficos, como vimos en el caso de la distribucion Normal.

>> subplot(1,2,1)

>> plot(x5,y5,., x10,y10, +, x20, y20, *);

>> legend(n=5, p=0.2, n=10, p=0.2, n=20, p=0.2)

>> subplot(1,2,2)

>> plot(x100, y1, ., x100, y2, +, x100, y3, *);

>> legend(n = 100, p = 0.1, n = 100, p = 0.5, n = 100, p = 0.8);

0 5 10 15 200

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

n=5, p=0.2n=10, p=0.2n=20, p=0.2

0 20 40 60 80 1000

0.02

0.04

0.06

0.08

0.1

0.12

0.14

n = 100 , p = 0.1n = 100 , p = 0.5n = 100 , p = 0.8

4. Sabiendo que la distribucion Bin(n, p) es la suma de n variables aleatorias Bernouilli inde-


pendientes:Binom(n, p) = Bern(p) + ...+ Bern(p)

n veces

Crea una funcion en MATLAB/Octave para generar numeros aleatorios de la distribucionbinomial.

%% creamos la funcion na_binomial.m

function y= na_binomial(n_dat, n, p) % n_dat es el

% no de datos binomiales % que se van a generar

prob = rand(n,n_dat); % realizamos n experimentos

exitos = prob < p; % si se cumple, entonces se ha producido un exito

y =sum(exitos); % sumamos todos los exitos producidos

%%

En el Command Window:

% Obtenemos 10 numeros binomiales obtenidos de sumar el no

% de exitos de n intentos con una probabilidad de exito p

>> y = na_binomial(10,n,p)

Numericamente, podemos verificar las propiedades de la media y la varianza de la distri-bucion binomial:

E[X] = n pVar[X] = n p q

>> n=100; p=0.1;

>> y = na_binomial(10,n,p)

>> mean(y) % aprox. 10

>> var(y) % aprox. 9

3.1. Generacion de variables aleatorias

Como se indico en la practica anterior, el metodo de la inversa nos permite generar variablesaleatorias contnuas a partir de la funcion de distribucion siempre que esta admita inversa.


Caso Weibull: Sea X una variable aleatoria Weibull2, X Weibull(, ), con funcion dedistribucion:

FX(x) =

{0 x > 0

1 e(x/) 0 xdonde y son parametros dados. Indica como generaras en MATLAB/Octave variablesaleatorias Weibull.

X Mediante el metodo de la inversa, tenemos que considerando u = FX(x),

1 e(x/) = ue(x/)

= 1 u

(x/) = ln(1 u)x/ = [ ln(1 u)] 1x = [ ln(1 u)] 1

Por tanto:

>> u = rand(n,1);

>> x = beta*(-log(1-u))^(1/alpha);

4. Funcion Q

La funcion Q se define como la complementaria a la funcion de distribucion de la N (0, 1), esdecir como

Q(x) = P (X > x) siendo X N (0, 1)Supongamos una funcion Q tiene por cotas superiores (cs1) y (cs2) y cota inferior (ci) dadas por

Q(x) 12e

x2

2 para todo x 0 (cs1)

Q(x) 0 (cs2)

Q(x) >12pix

(1 1

x2

)e

x2

2 para todo x > 1 (ci)

2NOTA: Esta distribucion se aplica en los analisis de fiabilidad de sistemas, para establecer, por ejemplo el tiempode vida de un componente hasta que se produce un fallo.


a) La funcion de MATLAB normcdf(x,MU,SIGMA) devuelve la funcion de distribucion acumu-lada de una N (, ). Como representaras en MATLAB la funcion Q?

>> Q = 1 - normcdf(x,0,1)

b) Escribe el codigo en MATLAB para representar graficamente la funcion Q en el intervalo(0, 5] junto con las cotas superiores cs1 y cs2 y la cota inferior ci. El resultado que se ha deobtener sera similar al de la figura 1.

NOTA:

Representa x (0, 5] como x=[0.01:0.01:5].Para la cota inferior ci, observa que existe una asntota vertical en x = 1, por tanto,define xi=[1.01:0.01:5] tan solo para esta cota.

Utiliza una escala logartmica para representar graficamente las cotas.

>>x1= [1.01:0.01:5];

>>for i=1:length(x1)

ci(i) = (1/x1(i)*sqrt(2*pi))*(1- (1/x1(i)^2))*exp(-0.5*x1(i)^2);

end

>>x = [0.01:0.01:5];

>>for i=1:length(x)

cs1(i)=0.5*exp(-0.5*x(i)^2);

cs2(i)=(1/(x(i)*sqrt(2*pi))) * exp(-0.5*x(i)^2);

end

>>q = 1- normcdf(x,0,1);

>> plot(x, log(q), black)

>> hold on

>> plot(x, log(cs1), r)

>> plot(x, log(cs2), g)

>> plot(x1, log(ci), m)

>> hold off


0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-20

-15

-10

-5

0

5

10log(q)log(cs1)log(cs2)log(ci)

Figura 1: Representacion grafica de la funcion Q con las cotas en escala logartmica log(cs1),log(cs2) y log(ci). La lnea vertical representa la asntota en x = 1.

5. Ejercicios de interes

1. Si se sabe que una v.a Y del tipo Cauchy se puede obtener como Y = tg (X) con X unaU(pi2 , pi2 ), como se puede aprovechar esta informacion para generar una Cauchy?Si se sabe dicha informacion no hace falta recurrir al metodo de la inversa de la funcion

de distribucion. Simplemente podemos:

>> x=unifrnd(-pi/2,pi/2, 100,1); % generar uniformes en intervalo

% (-pi/2,pi/2)

>> y=tan(x); % aplicar la tangente

2. Sea X N (0, = 3). Se construye un crculo con radio un valor |x| de la v.a. anterior.Calcular en terminos de la funcion Q la probabilidad de que el crculo generado tenga un areamayor o igual que pi. Evaluar dicha probabilidad con MATLAB/ Octave de forma exacta ypor simulacion.


Sea C la v.a. Area del crculo, se trata de calcular

P (C pi) = P (piX2 pi) = P (X2 1) = P (X 1) + P (X 1)= P

(X 0

3 1 0

3

)+ P

(X 0

3 1 0

3

)= 2Q

(1

3

)Para evaluarla con MATLAB/Octave de forma exacta:

>> x=2*(1-normcdf(1/3))

0.7389

Para resolverlo con MATLAB/Octave de forma simulada:

>> x=normrnd(0,3,1000,1);

>> area=pi*x.>2;

>> c=(area>=pi);

>> prob=sum(c)/1000

0.7080

3. Sea X la v.a. horas que se dedica a realizar una actividad, cuya funcion de densidad es

f (x) =

{14 (x+ 1) 0 < x < 2

0 resto

Se pide:

a) Calcular, de forma analtica y con MATLAB/Octave, la probabilidad de que el tiempoempleado sea superior a una hora y media.

b) Si se realizan 10 actividades segun la v.a.X, calcular, de forma analtica y con MATLAB/Octave,la probabilidad de que exactamente en tres de ellas, el tiempo que se emplee en realizarcada una sea superior a una hora y media.


a) De forma analtica, se tiene que:

P (X > 1,5) =

+x=1,5

f (x) dx =

2x=1,5

1

4(x+ 1) dx =

1

4

1

2

[x2]21,5

+1

4[x]21,5

=1

8

(4 9

4

)+

1

4

(2 3

2

)=

7

32+

1

8=

11

32= 0,3438

Para resolverlo con MATLAB/Octave, hay que utilizar el metodo de la inversa, con loque hay que calcular previamente la funcion de distribucion, que es

F (x) =

0 x < 0 x

014 (y + 1) dy =

18x

2 + 14x 0 x < 21 2 x

Y ahora hay que invertir la funcion de distribucion:

1

8x2 +

1

4x = u

x2 + 2x 8u = 0

x =24 + 32u

2

Dada la definicion de la v.a. X la solucion negativa no es valida, con lo que el proce-dimiento para generarla es:

u U (0, 1)

x =2 +4 + 32u

2

El codigo MATLAB/Octave para resolver este apartado es:

>> u=rand(1000,1);

>> x=(-2+sqrt(4+32*u))/2;

>> prob=sum(x>1.5)/1000

0.33800

b) De forma analtica, si se define Y como la v.a. numero de actividades entre las 10realizadas en las que el tiempo que se emplea es superior a una hora y media, se tieneque Y Bin (n = 10, p = 1132) hay que calcular

P (Y = 3) =

(10

3

)(11

32

)3(2132

)7= 0, 2555


El codigo MATLAB/Octave para resolver este apartado es:

>> u=rand(1000,10);

>> x=(-2+sqrt(4+32*u))/2;

>> conta=(x>1.5);

>> suma=sum(conta,2);

>> prob=sum(suma==3)/1000

0.24500

4. Ex. FEB 2004 Ing. Tel. (C1a) (link). El objetivo de este problema es analizar un canal decomunicaciones. Cuando el canal transmite un 1, el receptor recibe una variable aleatoria quesigue una distribucion Normal de media 1 y varianza 0, 5. Si el canal binario transmite un2, el receptor recibe una variable aleatoria Normal con media 2 y varianza 0, 5. Sea P (1) laprobabilidad de transmitir un 1.

a) Si P (1) = 0, 75. Cual es la probabilidad de que un 1 haya sido transmitido cuando elreceptor ha recibido una senal superior a 2?

T1 = transmitir un 1, P (T1) = 0,75 , R|T1 N(1, 0,5)T2 = transmitir un 2, P (T2) = 0,25 , R|T2 N(2, 0,5)

P (T1|R > 2) = P (R > 2|T1)P (T1)P (R > 2)

=

=P (R > 2|T1)P (T1)

P (R > 2|T1)P (T1) + P (R > 2|T2)P (T2)

Sean Z1 =R|T11

0,5y Z2 =

R|T220,5

,

P (T1|R > 2) = P (Z1 > 1,41) 0,75P (Z1 > 1,41) 0,75 + P (Z2 > 0) 0,25

=0,0793 0,75

0,0793 0,75 + 0,5 0,25 = 0,3224

En MATLAB/Octave, podemos aproximar la probabilidad mediante:

>> n=10000;

>> u=rand(n,1); % simulamos n experimientos aleatorios

>> p=0.75; % es la probabilidad de transmitir un 1

% transmision del mensaje por el canal

>> t=1*(up);

% recepcion del mensaje

>> r=(t==1).*normrnd(1,sqrt(0.5),n,1)+(t==2).*normrnd(2,sqrt(0.5),n,1);

>> r=2*(r>2)+1*(r

% creamos una tabla de doble entrada transmitidos/recibidos (SOLO EN MATLAB)

>> a=crosstab(t,r)

>> a =

6894 599 % a(1,1) a(1,2)

1258 1249 % a(2,1) a(2,2)

% nos piden la probabilidad P(t==1|r>2), en la simulacion

>> a(1,2)/(a(1,2)+a(2,2))

0.3241

5. Ex. SEP 2007 Ing. Tec. Tel. (P2b). (link). La duracion en das de un tipo de sensores sigue unmodelo Weibull con

F (t) = 1 exp((t

)1/2)

f (t) =1/2

1/2t1/2 exp

((t

)1/2)

con > 0.

a) Se supone que se tiene una caja de 60 sensores sin usar cuya duracion siguen el modeloWeibull con = 14 das que verifica que E [T ] =

12 das y V [T ] =

54 das

2. Si secomienza con un sensor de dicha caja y se va reemplazando instantaneamente segun sevaya fundiendo con un sensor de la misma caja, cual es la probabilidad de que cuandohaya fallado el ultimo sensor de la caja hayan pasado menos de 47 das?

Utilizar alguno de los siguientes valores de la funcion

x 1,64 1,64 1,96 1,96 47 47Q (x) 0,05 0,95 0,025 0,975 ' 0 ' 1

b) Comentar lnea a lnea el siguiente codigo de MATLAB y decir el valor aproximado quetomara prob.

>> u=rand(60,1000);

>> w=0.25*(-log(1-u)).^(1/0.5);

>> s=sum(w);

>> prob=sum(s

a) Si se define la v.a. Ti como el tiempo hasta que falla el sensor i, se tiene que se pide

P

(60i=1

Ti 47)

Como las Ti verifican las condiciones del teorema central del lmite, se tiene que

60i=1

Ti N(

601

2,

60

5

4

)

Si se estandariza

P

(60i=1

Ti 47)

= P

60i=1 Ti 60126054

modelos probabilisticos

Documents

Transcript of modelos probabilisticos