Search

Text Color

Margin Size

Font Type

Enable Dyslexic Font

Inicio
Estadísticas
Teoría de Probabilidad
Probabilidad Aplicada (Pfeiffer)
17: Apéndices
17.1: Apéndice A a Probabilidad Aplicada- Directorio de m-funciones y m-procedimientos

17.1: Apéndice A a Probabilidad Aplicada- Directorio de m-funciones y m-procedimientos

Última actualización

31 oct 2022
Guardar como PDF
- 17: Apéndices
- 17.2: Apéndice B a Probabilidad Aplicada- algunas ayudas matemáticas

Paul Pfeiffer
Rice University

$\newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$

$\newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$

$\newcommand{\id}{\mathrm{id}}$ $\newcommand{\Span}{\mathrm{span}}$

( \newcommand{\kernel}{\mathrm{null}\,}\) $\newcommand{\range}{\mathrm{range}\,}$

$\newcommand{\RealPart}{\mathrm{Re}}$ $\newcommand{\ImaginaryPart}{\mathrm{Im}}$

$\newcommand{\Argument}{\mathrm{Arg}}$ $\newcommand{\norm}[1]{\| #1 \|}$

$\newcommand{\inner}[2]{\langle #1, #2 \rangle}$

$\newcommand{\Span}{\mathrm{span}}$

$\newcommand{\id}{\mathrm{id}}$

$\newcommand{\Span}{\mathrm{span}}$

$\newcommand{\kernel}{\mathrm{null}\,}$

$\newcommand{\range}{\mathrm{range}\,}$

$\newcommand{\RealPart}{\mathrm{Re}}$

$\newcommand{\ImaginaryPart}{\mathrm{Im}}$

$\newcommand{\Argument}{\mathrm{Arg}}$

$\newcommand{\norm}[1]{\| #1 \|}$

$\newcommand{\inner}[2]{\langle #1, #2 \rangle}$

$\newcommand{\Span}{\mathrm{span}}$ $\newcommand{\AA}{\unicode[.8,0]{x212B}}$

$\newcommand{\vectorA}[1]{\vec{#1}} % arrow$

$\newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow$

$\newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$

$\newcommand{\vectorC}[1]{\textbf{#1}}$

$\newcommand{\vectorD}[1]{\overrightarrow{#1}}$

$\newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}}$

$\newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}}$

$\newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$

$\newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$

$\newcommand{\avec}{\mathbf a}$

$\newcommand{\bvec}{\mathbf b}$

$\newcommand{\cvec}{\mathbf c}$

$\newcommand{\dvec}{\mathbf d}$

$\newcommand{\dtil}{\widetilde{\mathbf d}}$

$\newcommand{\evec}{\mathbf e}$

$\newcommand{\fvec}{\mathbf f}$

$\newcommand{\nvec}{\mathbf n}$

$\newcommand{\pvec}{\mathbf p}$

$\newcommand{\qvec}{\mathbf q}$

$\newcommand{\svec}{\mathbf s}$

$\newcommand{\tvec}{\mathbf t}$

$\newcommand{\uvec}{\mathbf u}$

$\newcommand{\vvec}{\mathbf v}$

$\newcommand{\wvec}{\mathbf w}$

$\newcommand{\xvec}{\mathbf x}$

$\newcommand{\yvec}{\mathbf y}$

$\newcommand{\zvec}{\mathbf z}$

$\newcommand{\rvec}{\mathbf r}$

$\newcommand{\mvec}{\mathbf m}$

$\newcommand{\zerovec}{\mathbf 0}$

$\newcommand{\onevec}{\mathbf 1}$

$\newcommand{\real}{\mathbb R}$

$\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}$

$\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}$

$\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}$

$\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}$

$\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}$

$\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}$

$\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}$

$\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}$

$\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}$

$\newcommand{\laspan}[1]{\text{Span}\{#1\}}$

$\newcommand{\bcal}{\cal B}$

$\newcommand{\ccal}{\cal C}$

$\newcommand{\scal}{\cal S}$

$\newcommand{\wcal}{\cal W}$

$\newcommand{\ecal}{\cal E}$

$\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}$

$\newcommand{\gray}[1]{\color{gray}{#1}}$

$\newcommand{\lgray}[1]{\color{lightgray}{#1}}$

$\newcommand{\rank}{\operatorname{rank}}$

$\newcommand{\row}{\text{Row}}$

$\newcommand{\col}{\text{Col}}$

$\renewcommand{\row}{\text{Row}}$

$\newcommand{\nul}{\text{Nul}}$

$\newcommand{\var}{\text{Var}}$

$\newcommand{\corr}{\text{corr}}$

$\newcommand{\len}[1]{\left|#1\right|}$

$\newcommand{\bbar}{\overline{\bvec}}$

$\newcommand{\bhat}{\widehat{\bvec}}$

$\newcommand{\bperp}{\bvec^\perp}$

$\newcommand{\xhat}{\widehat{\xvec}}$

$\newcommand{\vhat}{\widehat{\vvec}}$

$\newcommand{\uhat}{\widehat{\uvec}}$

$\newcommand{\what}{\widehat{\wvec}}$

$\newcommand{\Sighat}{\widehat{\Sigma}}$

$\newcommand{\lt}{<}$

$\newcommand{\gt}{>}$

$\newcommand{\amp}{&}$

$\definecolor{fillinmathshade}{gray}{0.9}$

Usamos el término m-function para designar una función definida por el usuario como distinta de las funciones básicas de MATLAB que forman parte del paquete MATLAB. Por ejemplo, la función m minterm produce el vector minterm especificado. Un procedimiento m (o a veces un procedimiento) es un archivo m que contiene un conjunto de comandos MATLAB que llevan a cabo un conjunto prescrito de operaciones. Generalmente, estos solicitarán (o asumirán) ciertos datos sobre los cuales se lleva a cabo el trámite. Utilizamos el término m-programa para referirnos a una función m o a un procedimiento m.

Además de los programas m hay una colección de archivos m con datos correctamente formateados que se pueden ingresar al espacio de trabajo llamando al archivo.

Aunque los programas m fueron escritos para MATLAB versión 4.2, funcionan para las versiones 5.1, 5.2 y 7.04. Estas últimas versiones ofrecen algunas características nuevas que pueden hacer más eficiente la implementación de algunos de los programas m, y que hacen posibles algunos nuevos. Con una excepción (así señalada), estos no se exploran en esta colección.

Características de MATLAB

La utilización de los recursos de MATLAB es posible gracias a un análisis sistemático de algunas características del modelo básico de probabilidad. En particular, se explota el análisis minterm de combinaciones lógicas (o booleanas) de eventos y el análisis de la estructura de variables aleatorias simples con la ayuda de funciones indicadoras y análisis minterm.

Una serie de características estándar de MATLAB se utilizan ampliamente. Además del álgebra matricial estándar, utilizamos:

Aritmética de matriz. Esto implica cálculos elemento por elemento. Por ejemplo, si a, b son matrices del mismo tamaño, entonces a.*b es la matriz obtenida multiplicando los elementos correspondientes en las dos matrices para obtener una nueva matriz del mismo tamaño.
Operaciones relacionales, como menor que, igual, etc. para obtener matrices cero-uno con unas en posiciones de elementos donde se cumplen las condiciones.
Operaciones lógicas en matrices cero-uno utilizando operadores lógicos y, o, y no, así como ciertas funciones relacionadas como cualquiera, todas, not, find, etc. Nota. Las operaciones relacionales y las operaciones lógicas producen matrices cero-uno, llamadas matrices lógicas, que MATLAB trata de manera diferente a las matrices numéricas cero-uno. Una matriz rectangular en la que algunas filas son matrices lógicas pero otras no se trata como una matriz numérica. Cualquier matriz rectangular cero-uno se puede convertir en una matriz numérica (matriz) mediante el comando A = unos (tamaño (A)). *A,
Ciertas funciones de MATLAB, como meshgrid, sum, cumsum, prod, cumprod se utilizan repetidamente. La función punto para el producto punto no funciona si cualquiera de las matrices es una matriz lógica. Si uno de los pares es numérico, funcionará el comando C = A*B'.

Bloques de construcción auxiliares definidos por el usuario

csort.m

Descripción del Código:

Una de las más útiles es una operación especial de clasificación y consolidación implementada en la función m csort. Un problema estándar surge cuando cada uno de un conjunto no distinto de valores tiene una probabilidad asociada. Para obtener la distribución, es necesario ordenar los valores y sumar las probabilidades asociadas a cada valor distinto. La siguiente función m logra estas operaciones: function [t, p] = csort (T, P). T y P son matrices con el mismo número de elementos. Se ordenan los valores de T y se consolidan valores idénticos; se suman valores de P correspondientes a valores idénticos de T. Una serie de funciones y procedimientos derivados utilizan csort. Los dos siguientes son útiles.

Contestar

function [t,p] = csort(T,P)
% CSORT  [t,p] = csort(T,P) Sorts T, consolidates P
% Version of 4/6/97
% Modified to work with Versions 4.2 and 5.1, 5.2
% T and P matrices with the same number of elements
% The vector T(:)' is sorted:
%   * Identical values in T are consolidated;
%   * Corresponding values in P are added.
T = T(:)';
n = length(T);
[TS,I] = sort(T);
d  = find([1,TS(2:n) - TS(1:n-1) >1e-13]); % Determines distinct values
t  = TS(d);                                % Selects the distinct values
m  = length(t) + 1;
P  = P(I);                                 % Arranges elements of P
F  = [0 cumsum(P(:)')];
Fd = F([d length(F)]);                     % Cumulative sums for distinct values
p  = Fd(2:m) - Fd(1:m-1);                  % Separates the sums for these values

distint.m

Descripción del Código:

función distint.m y = distinto (T) determina y ordena los distintos miembros de la matriz $T$ .

Contestar

function y = distinct(T)
% DISTINCT y = distinct(T) Disinct* members of T
% Version of 5/7/96  Rev 4/20/97 for version 4 & 5.1, 5.2
% Determines distinct members of matrix T.
% Members which differ by no more than 10^{-13}
% are considered identical.  y is a row
% vector of the distinct members.
TS = sort(T(:)');
n  = length(TS);
d  = [1  abs(TS(2:n) - TS(1:n-1)) >1e-13];
y  = TS(find(d));

freq.m

Descripción del Código:

freq.m ordena los distintos miembros de una matriz, cuenta el número de ocurrencias de cada valor y calcula las frecuencias relativas acumulativas.

Contestar

% FREQ file freq.m Frequencies of members of matrix
% Version of 5/7/96
% Sorts the distinct members of a matrix, counts
% the number of occurrences of each value, and
% calculates the cumulative relative frequencies.
T  = input('Enter matrix to be counted  ');
[m,n] = size(T);
[t,f] = csort(T,ones(m,n));
p = cumsum(f)/(m*n);
disp(['The number of entries is  ',num2str(m*n),])
disp(['The number of distinct entries is  ',num2str(length(t)),] )
disp(' ')
dis = [t;f;p]';
disp('    Values     Count   Cum Frac')
disp(dis)

dsum.m

Descripción del Código:

dsum.m función y = dsum (v, w) determina y ordena los distintos elementos entre las sumas de pares de elementos de los vectores de fila v y w.

Contestar

function y = dsum(v,w)
% DSUM y = dsum(v,w) Distinct pair sums of elements
% Version of 5/15/97
% y is a row vector of distinct
% values among pair sums of elements
% of matrices v, w.
% Uses m-function distinct
[a,b] = meshgrid(v,w);
t = a+b;
y = distinct(t(:)');

rep.m

Descripción del Código:

función rep.m y = rep (A, m, n) replica la matriz A, m veces verticalmente y n veces horizontalmente. Esencialmente lo mismo que la función repmat en MATLAB versión 5, publicada en diciembre de 1996.

Contestar

function y = rep(A,m,n)
% REP y = rep(A,m,n) Replicates matrix A
% Version of 4/21/96
% Replicates A,
% m times vertically,
% n times horizontally
% Essentially the same as repmat in version 5.1, 5.2
[r,c] = size(A);
R = [1:r]';
C = [1:c]';
v = R(:,ones(1,m));
w = C(:,ones(1,n));
y = A(v,w);

elrep.m

Descripción del Código:

elrep.m función y = elrep (A, m, n) replica cada elemento de A, $m$ veces verticalmente y $n$ veces horizontalmente.

Contestar

function y = elrep(A,m,n)
% ELREP y = elrep(A,m,n) Replicates elements of A
% Version of 4/21/96
% Replicates each element,
% m times vertically,
% n times horizontally
[r,c] = size(A);
R = 1:r;
C = 1:c;
v = R(ones(1,m),:);
w = C(ones(1,n),:);
y = A(v,w);

kronf.m

Descripción del Código:

kronf.m función y = kronf (A, B) determina el producto Kronecker de las matrices A, B Logra el mismo resultado para matrices completas que la función kron MATLAB.

Contestar

function y = kronf(A,B)
% KRONF y = kronf(A,B) Kronecker product
% Version of 4/21/96
% Calculates Kronecker product of full matrices.
% Uses m-functions elrep and rep
% Same result for full matrices as kron for version 5.1, 5.2
[r,c] = size(B);
[m,n] = size(A);
y = elrep(A,r,c).*rep(B,m,n);

colcopy.m

Descripción del Código:

colcopy.m función y = colcopy (v, n) trata el vector de fila o columna v como un vector de columna y hace una matriz con $n$ columnas de v.

Contestar

function y = colcopy(v,n)
% COLCOPY y = colcopy(v,n)  n columns of v
% Version of 6/8/95 (Arguments reversed 5/7/96)
% v a row or column vector
% Treats v as column vector
% and makes n copies
% Procedure based on "Tony's trick"
[r,c] = size(v);
if r == 1
  v = v';
end
y = v(:,ones(1,n));

colcopyi.m

Descripción del Código:

colcopyi.m función y = colcopyi (v, n) trata el vector de fila o columna v como un vector de columna, invierte el orden de los elementos y crea una matriz con n columnas del vector invertido.

Contestar

function y = colcopyi(v,n)
% COLCOPYI y = colcopyi(v,n) n columns in reverse order
% Version of 8/22/96
% v a row or column vector.
% Treats v as column vector,
% reverses the order of the
% elements, and makes n copies.
% Procedure based on "Tony's trick"
N = ones(1,n);
[r,c] = size(v);
if r == 1
  v = v(c:-1:1)';
else
  v = v(r:-1:1);
end
y = v(:,N);

rowcopy.m

Descripción del Código:

rowcopy.m función y = rowcopy (v, n) trata el vector de fila o columna v como un vector de fila y crea una matriz con $n$ filas de v.

Contestar

function y = rowcopy(v,n)
% ROWCOPY y = rowcopy(v,n)  n rows of v
% Version of 5/7/96
% v a row or column vector
% Treats v as row vector
% and makes n copies
% Procedure based on "Tony's trick"
[r,c] = size(v);
if c == 1
  v = v';
end
y = v(ones(1,n),:);

repseq.m

Descripción del Código:

repseq.m función y = repseq (V, n) replica los $V$ $n$ tiempos vectoriales, horizontalmente si $V$ es un vector de fila y verticalmente si $V$ es un vector de columna.

Contestar

function y = repseq(V,n);
% REPSEQ y = repseq(V,n) Replicates vector V n times
% Version of 3/27/97
% n replications of vector V
% Horizontally if V a row vector
% Vertically if V a column vector
m = length(V);
s = rem(0:n*m-1,m)+1;
y = V(s);

total.m

Descripción del Código:

total.m Total de todos los elementos en una matriz, calculado por: total (x) = suma (suma (x)).

Contestar

function y = total(x)
% TOTAL y = total(x)
% Version of 8/1/93
% Total of all elements in matrix x.
y = sum(sum(x));

dispv.m

Descripción del Código:

dispv.m Las matrices $A, B$ se transponen y se muestran una al lado de la otra.

Contestar

function y = dispv(A,B)
% DISPV y = dispv(A,B) Transpose of A, B side by side
% Version of 5/3/96
% A, B are matrices of the same size
% They are transposed and displayed
% side by side.
y  = [A;B]';

roundn.m

Descripción del Código:

roundn.m función y = roundn (A, n) redondea la matriz A a $n$ decimales.

Contestar

function y = roundn(A,n);
% ROUNDN y = roundn(A,n)
% Version of 7/28/97
% Rounds matrix A to n decimals
y = round(A*10^n)/10^n;

arrep.m

Descripción del Código:

arrep.m función y = arrep (n, k) forma todos los arreglos, con repetición, de $k$ elementos de la secuencia $1: n$ .

Contestar

function y = arrep(n,k);
% ARREP  y = arrep(n,k);
% Version of 7/28/97
% Computes all arrangements of k elements of 1:n,
% with repetition allowed. k may be greater than n.
% If only one input argument n, then k = n.
% To get arrangements of column vector V, use
% V(arrep(length(V),k)).
N = 1:n;
if nargin == 1
  k = n;
end
y = zeros(k,n^k);
for i = 1:k
  y(i,:) = rep(elrep(N,1,n^(k-i)),1,n^(i-1));
end

Vectores y probabilidades minterm

El análisis de combinaciones lógicas de eventos (como conjuntos) se sistematiza mediante el uso de la expansión minterm. Esto conduce naturalmente a la noción de vectores minterm. Estos son vectores cero-uno que se pueden combinar mediante operaciones lógicas. La producción de los patrones básicos de minterm es esencial para una serie de operaciones. Los siguientes m-programas son elementos clave de varios otros programas.

minterm.m

Descripción del Código:

minterm.m función y = minterm (n, k) genera el $k$ th vector minterm en una clase de $n$ .

Contestar

function y = minterm(n,k)
% MINTERM y = minterm(n,k) kth minterm of class of n
% Version of 5/5/96
% Generates the kth minterm vector in a class of n
% Uses m-function rep
y = rep([zeros(1,2^(n-k)) ones(1,2^(n-k))],1,2^(k-1));

mintable.m

Descripción del Código:

mintable.m función y = mintable (n) genera una tabla de vectores minterm mediante el uso repetido de la función m minterm.

Contestar

function y = mintable(n)
% MINTABLE y = mintable(n)  Table of minterms vectors
% Version of 3/2/93
% Generates a table of minterm vectors
% Uses the m-function minterm
y = zeros(n,2^n);
for i = 1:n
    y(i,:) = minterm(n,i);
end

minvec3.m

Descripción del Código:

minvec3.m establece vectores minterm básicos A, B, C, A $^c$ , B $^c$ , C $^c$ , para la clase $\{A, B, C\}$ . (De manera similar para minvec4.m, minvec5.m, etc.)

Contestar

% MINVEC3  file minvec3.m Basic minterm vectors
% Version of 1/31/95
A = minterm(3,1);
B = minterm(3,2);
C = minterm(3,3);
Ac = ~A;
Bc = ~B;
Cc = ~C;
disp('Variables are A, B, C, Ac, Bc, Cc')
disp('They may be renamed, if desired.')

minmap

Descripción del Código:

minmap función y = minmap (pm) reforma un vector de fila o columna pm de probabilidades minterm en formato de mapa minterm.

Contestar

function y = minmap(pm)
% MINMAP y = minmap(pm) Reshapes vector of minterm probabilities 
% Version of 12/9/93 
% Reshapes a row or column vector pm of minterm 
% probabilities into minterm map format
m = length(pm);
n = round(log(m)/log(2));
a = fix(n/2);
if m ~= 2^n
  disp('The number of minterms is incorrect')
else
  y = reshape(pm,2^a,2^(n-a));
end

binario.m

Descripción del Código:

binary.m función y = binario (d, n) convierte una matriz d de enteros no negativos de punto flotante en una matriz de equivalentes binarios, uno en cada fila. Adaptado de m-funciones escritas por Hans Olsson y por Simon Cooke. Cada fila de matriz se puede convertir en una cadena no espaciada de ceros y unos por el dispositivo ys = setstr (y + '0').

Contestar

function y = binary(d,n)
% BINARY y = binary(d,n) Integers to binary equivalents
% Version of 7/14/95 
% Converts a matrix d of floating point, nonnegative
% integers to a matrix of binary equivalents. Each row
% is the binary equivalent (n places) of one number.
% Adapted from the programs dec2bin.m, which shared
% first prize in an April 95 Mathworks contest. 
% Winning authors: Hans Olsson from Lund, Sweden,
% and Simon Cooke from Glasgow, UK.
% Each matrix row may be converted to an unspaced string 
% of zeros and ones by the device:  ys = setstr(y + '0').
if nargin < 2,  n = 1; end     % Allows omission of argument n
[f,e] = log2(d); 
n = max(max(max(e)),n);      
y = rem(floor(d(:)*pow2(1-n:0)),2);

mincalc.m

Descripción del Código:

mincalc.m El procedimiento m mincalc determina las probabilidades minterm a partir de datos adecuados. Para una discusión sobre el formato de los datos y ciertos problemas, ver 2.6.

Contestar

% MINCALC file mincalc.m Determines minterm probabilities
% Version of 1/22/94 Updated for version 5.1 on  6/6/97
% Assumes a data file which includes
%  1. Call for minvecq to set q basic minterm vectors, each (1 x 2^q)
%  2. Data vectors DV = matrix of md data Boolean combinations of basic sets--
%     Matlab produces md minterm vectors-- one on each row.
%     The first combination is always A|Ac (the whole space)
%  3. DP = row matrix of md data probabilities.
%     The first probability is always 1.
%  4. Target vectors TV = matrix of mt target Boolean combinations.
%     Matlab produces a row minterm vector for each target combination.
%     If there are no target combinations, set TV = [];
[md,nd] = size(DV);
ND = 0:nd-1;
ID = eye(nd);                 % Row i is minterm vector i-1
[mt,nt] = size(TV);
MT = 1:mt;
rd = rank(DV);
if rd < md                    
   disp('Data vectors are NOT linearly independent')
  else
   disp('Data vectors are linearly independent')
end
% Identification of which minterm probabilities can be determined from the data
% (i.e., which minterm vectors are not linearly independent of data vectors)
AM = zeros(1,nd);
for i = 1:nd
  AM(i) = rd == rank([DV;ID(i,:)]);  % Checks for linear dependence of each
  end   
am = find(AM);                             % minterm vector
CAM = ID(am,:)/DV;     % Determination of coefficients for the available minterms
pma = DP*CAM';                % Calculation of probabilities of available minterms
PMA = [ND(am);pma]';
if sum(pma < -0.001) > 0      % Check for data consistency
   disp('Data probabilities are INCONSISTENT')
else
% Identification of which target probabilities are computable from the data
CT = zeros(1,mt);
for j = 1:mt
  CT(j) = rd == rank([DV;TV(j,:)]);
  end
ct  = find(CT);
CCT = TV(ct,:)/DV;            % Determination of coefficients for computable targets
ctp = DP*CCT';                % Determination of probabilities
disp(' Computable target probabilities')
disp([MT(ct); ctp]')
end                           % end for "if sum(pma < -0.001) > 0"
disp(['The number of minterms is ',num2str(nd),])
disp(['The number of available minterms is ',num2str(length(pma)),])
disp('Available minterm probabilities are in vector pma')
disp('To view available minterm probabilities, call for PMA')

mincalct.m

Descripción del Código:

mincalct.m Modificación de mincalc. Asume que mincalc se ha ejecutado, pide nuevos vectores objetivo y realiza los mismos cálculos que mincalc.

Contestar

% MINCALCT file mincalct.m  Aditional target probabilities
% Version of 9/1/93  Updated for version 5 on 6/6/97
% Assumes a data file which includes
%  1. Call for minvecq to set q basic minterm vectors.
%  2. Data vectors DV.  The first combination is always A|Ac.
%  3. Row matrix DP of data probabilities. The first entry is always 1.
TV = input('Enter matrix of target Boolean combinations  ');
[md,nd] = size(DV);
[mt,nt] = size(TV);
MT = 1:mt;
rd = rank(DV);
CT = zeros(1,mt);   % Identification of computable target probabilities
for j = 1:mt
  CT(j) = rd == rank([DV;TV(j,:)]);
end
ct  = find(CT);
CCT = TV(ct,:)/DV;  % Determination of coefficients for computable targets
ctp = DP*CCT';      % Determination of probabilities
disp(' Computable target probabilities')
disp([MT(ct); ctp]')

Eventos independientes

minprob.m

Descripción del Código:

minprob.m función y = minprob (p) calcula las probabilidades minterm para las probabilidades básicas en fila o columna vector p. Utiliza las funciones m mintable, colcopy.

Contestar

function y = minprob(p)
% MINPROB y = minprob(p) Minterm probs for independent events
% Version of 4/7/96
% p is a vector [P(A1) P(A2) ... P(An)], with
% {A1,A2, ... An} independent.
% y is the row vector of minterm probabilities
% Uses the m-functions mintable, colcopy
n = length(p);
M = mintable(n);
a = colcopy(p,2^n);          % 2^n columns, each the vector p
m = a.*M + (1 - a).*(1 - M); % Puts probabilities into the minterm 
                             % pattern on its side (n by 2^n)
y = prod(m);                 % Product of each column of m

imintest.m

Descripción del Código:

imintest.m función y = imintest (pm) comprueba probabilidades minterm de independencia.

Contestar

function y = imintest(pm)
% IMINTEST y = imintest(pm) Checks minterm probs for independence
% Version of 1/25//96
% Checks minterm probabilities for independence
% Uses the m-functions mintable and minprob
m = length(pm);
n = round(log(m)/log(2));
if m ~= 2^n
  y = 'The number of minterm probabilities is incorrect';
else
P = mintable(n)*pm';
pt = minprob(P');
a = fix(n/2);
s = abs(pm - pt) > 1e-7;
  if sum(s) > 0
    disp('The class is NOT independent')
    disp('Minterms for which the product rule fails')
    y = reshape(s,2^a,2^(n-a));
  else
    y = 'The class is independent';
  end
end

ikn.m

Descripción del Código:

ikn.m función y = ikn (P, k) determina la probabilidad de la ocurrencia de exactamente $k$ de los eventos $n$ independientes cuyas probabilidades están en fila o columna vector P
(k puede ser un vector de fila o columna de enteros no negativos menor o igual a $n$ ).

Contestar

function y = ikn(P,k)
% IKN y = ikn(P,k) Individual probabilities of k of n successes
% Version of 5/15/95
% Uses the m-functions mintable, minprob, csort
n = length(P);
T = sum(mintable(n));  % The number of successes in each minterm
pm = minprob(P);       % The probability of each minterm
[t,p] = csort(T,pm);   % Sorts and consolidates success numbers
                       % and adds corresponding probabilities
y = p(k+1);

ckn.m

Descripción del Código:

ckn.m función y = ckn (P, k) determina la probabilidad de la ocurrencia de $k$ o más de los eventos $n$ independientes cuyas probabilidades están en el vector de fila o columna P ( $k$ puede ser un vector de fila o columna)

Contestar

function y = ckn(P,k)
% CKN y = ckn(P,k) Probability of k or more successes
% Version of 5/15/95
% Probabilities of k or more of n independent events
% Uses the m-functions mintable, minprob, csort
n = length(P);
m = length(k);
T = sum(mintable(n));  % The number of successes in each minterm
pm = minprob(P);       % The probability of each minterm
[t,p] = csort(T,pm);   % Sorts and consolidates success numbers
                       % and adds corresponding probabilities
for i = 1:m            % Sums probabilities for each k value
  y(i) = sum(p(k(i)+1:n+1));
end

paralelo.m

Descripción del Código:

parallel.m función y = paralelo (p) determina la probabilidad de una combinación paralela de los eventos independientes cuyas probabilidades están en fila o columna vector p.

Contestar

function y = parallel(p)
% PARALLEL y = parallel(p) Probaaability of parallel combination
% Version of 3/3/93
% Probability of parallel combination.
% Individual probabilities in row matrix p.
y = 1 - prod(1 - p);

Probabilidad condicional e idependencia condicional

bayes.m

Descripción del Código:

bayes.m produce una inversión bayesiana de probabilidades condicionales. La entrada consiste en $P(E|A_i)$ y $P(A_i)$ para una clase disjunta $\{A_i: 1 \le i \le n\}$ cuya unión contiene $E$ . El procedimiento calcula $P(A_i|E)$ y $P(A_i|E^c)$ para $1 \le i \le n$ .

Contestar

% BAYES file bayes.m Bayesian reversal of conditional probabilities
% Version of 7/6/93 
% Input P(E|Ai) and P(Ai)
% Calculates P(Ai|E) and P(Ai|Ec)
disp('Requires input PEA = [P(E|A1) P(E|A2) ... P(E|An)]')
disp(' and PA = [P(A1) P(A2) ... P(An)]')
disp('Determines PAE  = [P(A1|E) P(A2|E) ... P(An|E)]')
disp('       and PAEc = [P(A1|Ec) P(A2|Ec) ... P(An|Ec)]')
PEA  = input('Enter matrix PEA of conditional probabilities  ');
PA   = input('Enter matrix  PA of probabilities  ');
PE   = PEA*PA';
PAE  = (PEA.*PA)/PE;
PAEc = ((1 - PEA).*PA)/(1 -  PE);
disp(' ')
disp(['P(E) = ',num2str(PE),])
disp(' ')
disp('    P(E|Ai)   P(Ai)     P(Ai|E)   P(Ai|Ec)')
disp([PEA; PA; PAE; PAEc]')
disp('Various quantities are in the matrices PEA, PA, PAE, PAEc, named above')

odds.m

Descripción del Código:

odds.m El procedimiento calcula las probabilidades posteriores para un perfil especificado $E$ . Asume que los datos han sido ingresados por el procedimiento oddsf u oddsp.

Contestar

% ODDS file odds.m  Posterior odds for profile
% Version of 12/4/93
% Calculates posterior odds for profile E
% Assumes data has been entered by oddsdf or oddsdp
E = input('Enter profile matrix E  ');
C =  diag(a(:,E))';       % aa = a(:,E) is an n by n matrix whose ith column
D =  diag(b(:,E))';       % is the E(i)th column of a.  The elements on the
                          % diagonal are b(i, E(i)), 1 <= i <= n
                          % Similarly for b(:,E)

R = prod(C./D)*(p1/p2);   % Calculates posterior odds for profile 
disp(' ')
disp(['Odds favoring Group 1:   ',num2str(R),])
if R > 1
  disp('Classify in Group 1')
else
  disp('Classify in Group 2')
end

oddsdf.m

Descripción del Código:

oddsdf.m Establece frecuencias de calibración para calcular las probabilidades posteriores.

Contestar

% ODDSDF file oddsdf.m Frequencies for calculating odds
% Version of 12/4/93
% Sets up calibrating frequencies 
% for calculating posterior odds
A = input('Enter matrix A of frequencies for calibration group 1  ');
B = input('Enter matrix B of frequencies for calibration group 2  ');
n = length(A(:,1));       % Number of questions (rows of A)
m = length(A(1,:));       % Number of answers to each question
p1 = sum(A(1,:));         % Number in calibration group 1
p2 = sum(B(1,:));         % Number in calibration group 2
a = A/p1;
b = B/p2;
disp(' ')                 % Blank line in presentation
disp(['Number of questions = ',num2str(n),])  % Size of profile
disp(['Answers per question = ',num2str(m),]) % Usually 3: yes, no, uncertain
disp(' Enter code for answers and call for procedure "odds"  ')
disp(' ')

oddsdp.m

Descripción del Código:

oddsdp.m Establece probabilidades condicionales para los cálculos de cuotas.

Contestar

  % ODDSDP file oddsdp.m Conditional probs for calculating posterior odds
% Version of 12/4/93
% Sets up conditional probabilities 
% for odds calculations
a  = input('Enter matrix A of conditional probabilities for Group 1  ');
b  = input('Enter matrix B of conditional probabilities for Group 2  ');
p1 = input('Probability p1 an individual is from Group 1  ');
n = length(a(:,1));
m = length(a(1,:));
p2 = 1 - p1;
disp(' ')                 % Blank line in presentation
disp(['Number of questions = ',num2str(n),])  % Size of profile
disp(['Answers per question = ',num2str(m),]) % Usually 3: yes, no, uncertain
disp(' Enter code for answers and call for procedure "odds"  ')
disp(' ')

Bernoulli y ensayos multinomiales

btdata.m

Descripción del Código:

btdata.m Establece el parámetro $p$ y el número $n$ de ensayos para generar secuencias de Bernoulli. Indicios para bt para generar los ensayos.

Contestar

  % BTDATA file btdata.m Parameters for Bernoulli trials
% Version of 11/28/92
% Sets parameters for generating Bernoulli trials
% Prompts for bt to generate the trials
n = input('Enter n, the number of trials  ');
p = input('Enter p, the probability of success on each trial  ');
disp(' ')
disp(' Call for bt')
disp(' ')

bt.m

Descripción del Código:

bt.m Genera secuencia de Bernoulli para parámetros establecidos por btdata. Calcula la frecuencia relativa de los “éxitos”.

Contestar

% BT file bt.m Generates Bernoulli sequence
% version of 8/11/95  Revised 7/31/97 for version 4.2 and 5.1, 5.2
% Generates Bernoulli sequence for parameters set by btdata
% Calculates relative frequency of 'successes'
clear SEQ;
B  = rand(n,1) <= p;         %  ones for random numbers <= p
F  = sum(B)/n;               % relative frequency of ones
N  = [1:n]';                 % display details
disp(['n = ',num2str(n),'   p = ',num2str(p),])
disp(['Relative frequency = ',num2str(F),])
SEQ = [N B];
clear N;
clear B;
disp('To view the sequence, call for SEQ')
disp(' ')

binomial.m

Descripción del Código:

binomial.m Utiliza ibinom y cbinom para generar tablas de las probabilidades binomiales individuales y acumulativas para parámetros especificados. Tenga en cuenta que para el cálculo en MATLAB suele ser mucho más conveniente y eficiente usar ibinom y/o cbinom.

Contestar

01% BINOMIAL file binomial.m Generates binomial tables
02% Version of 12/10/92  (Display modified 4/28/96)
03% Calculates a TABLE of binomial probabilities
04% for specified n, p, and row vector k,
05% Uses the m-functions ibinom and cbinom.
06n = input('Enter n, the number of trials ');
07p = input('Enter p, the probability of success ');
08k = input('Enter k, a row vector of success numbers ');
09y = ibinom(n,p,k);
10z = cbinom(n,p,k);
11disp(['    n = ',int2str(n),'    p = ' num2str(p)])
12H = ['    k         P(X = k)  P(X >= k)'];
13disp(H)
14disp([k;y;z]')

multinom.m

Descripción del Código:

multinom.m Distribución multinomial (pequeña $N, m$ ).

Contestar

% MULTINOM file multinom.m  Multinomial distribution
% Version of 8/24/96
% Multinomial distribution (small N, m)
N = input('Enter the number of trials  ');
m = input('Enter the number of types   ');
p = input('Enter the type probabilities  ');
M = 1:m;
T = zeros(m^N,N);
for i = 1:N
  a = rowcopy(M,m^(i-1));
  a = a(:);
  a = colcopy(a,m^(N-i));
  T(:,N-i+1) = a(:);        % All possible strings of the types
end
MT = zeros(m^N,m);
for i = 1:m
 MT(:,i) = sum(T'==i)';
end
clear T                     % To conserve memory
disp('String frequencies for type k are in column matrix MT(:,k)')
P = zeros(m^N,N);
for i = 1:N
  a = rowcopy(p,m^(i-1));
  a = a(:);
  a = colcopy(a,m^(N-i));
  P(:,N-i+1) = a(:);        % Strings of type probabilities
end
PS = prod(P');              % Probability of each string
clear P                     % To conserve memory
disp('String probabilities are in row matrix PS')

Algunos problemas de coincidencia

Cardmatch.m

Descripción del Código:

Cardmatch.m Muestreo para estimar la probabilidad de una o más coincidencias cuando se extrae una carta de cada una de las barajas $nd$ idénticas de $c$ cartas. Se especifica el número nsns de muestras.

Contestar

% CARDMATCH file cardmatch.m Prob of matches in cards from identical decks
% Version of 6/27/97
% Estimates the probability of one or more matches
% in drawing cards from nd decks of c cards each
% Produces a supersample of size n = nd*ns, where
% ns is the number of samples
% Each sample is sorted, and then tested for differences
% between adjacent elements.  Matches are indicated by
% zero differences between adjacent elements in sorted sample
c  = input('Enter the number c  of cards in a deck ');
nd = input('Enter the number nd of decks ');
ns = input('Enter the number ns of sample runs ');
X  = 1:c;                   % Population values
PX = (1/c)*ones(1,c);       % Population probabilities 
N  = nd*ns;                 % Length of supersample
U  = rand(1,N);             % Matrix of n random numbers
T  = dquant(X,PX,U);        % Supersample obtained with quantile function;
                            % the function dquant determines quantile
                            % function values of random number sequence U
ex = sum(T)/N;              % Sample average
EX = dot(X,PX);             % Population mean
vx = sum(T.^2)/N - ex^2;    % Sample variance
VX = dot(X.^2,PX) - EX^2;   % Population variance
A  = reshape(T,nd,ns);      % Chops supersample into ns samples of size nd
DS = diff(sort(A));         % Sorts each sample
m  = sum(DS==0)>0;          % Differences between elements in each sample
                            % Zero difference iff there is a match
pm = sum(m)/ns;             % Fraction of samples with one or more matches
Pm = 1 - comb(c,nd)*gamma(nd + 1)/c^(nd);  % Theoretical probability of match
disp('The sample is in column vector T')   % Displays of results
disp(['Sample average ex = ', num2str(ex),])
disp(['Population mean E(X) = ',num2str(EX),])
disp(['Sample variance vx = ',num2str(vx),])
disp(['Population variance V(X) = ',num2str(VX),])
disp(['Fraction of samples with one or more matches   pm = ', num2str(pm),])
disp(['Probability of one or more matches in a sample Pm = ', num2str(Pm),])

trialmatch.m

Descripción del Código:

trialmatch.m Estima la probabilidad de coincidencias en ensayos $n$ independientes a partir de distribuciones idénticas. El tamaño de la muestra y el número de ensayos deben mantenerse relativamente pequeños para evitar exceder la memoria disponible.

Contestar

% TRIALMATCH file trialmatch.m  Estimates probability of matches 
% in n independent trials from identical distributions
% Version of 8/20/97
% Estimates the probability of one or more matches 
% in a random selection from n identical distributions
% with a small number of possible values
% Produces a supersample of size N = n*ns, where
% ns is the number of samples.  Samples are separated.
% Each sample is sorted, and then tested for differences
% between adjacent elements.  Matches are indicated by
% zero differences between adjacent elements in sorted sample.
X  = input('Enter the VALUES in the distribution ');
PX = input('Enter the PROBABILITIES  ');
c  = length(X);
n  = input('Enter the SAMPLE SIZE n ');
ns = input('Enter the number ns of sample runs ');
N  = n*ns;                  % Length of supersample
U  = rand(1,N);             % Vector of N random numbers
T  = dquant(X,PX,U);        % Supersample obtained with quantile function;
                            %   the function dquant determines quantile
                            %   function values for random number sequence U
ex = sum(T)/N;              % Sample average
EX = dot(X,PX);             % Population mean
vx = sum(T.^2)/N - ex^2;    % Sample variance
VX = dot(X.^2,PX) - EX^2;   % Population variance
A  = reshape(T,n,ns);       % Chops supersample into ns samples of size n
DS = diff(sort(A));         % Sorts each sample
m  = sum(DS==0)>0;          % Differences between elements in each sample
                            % -- Zero difference iff there is a match
pm = sum(m)/ns;             % Fraction of samples with one or more matches
d  = arrep(c,n);
p  = PX(d);
p  = reshape(p,size(d));    % This step not needed in version 5.1
ds = diff(sort(d))==0;
mm = sum(ds)>0;
m0 = find(1-mm);
pm0 = p(:,m0);              % Probabilities for arrangements with no matches
P0 = sum(prod(pm0));      
disp('The sample is in column vector T')   % Displays of results
disp(['Sample average ex = ', num2str(ex),])
disp(['Population mean E(X) = ',num2str(EX),])
disp(['Sample variance vx = ',num2str(vx),])
disp(['Population variance V(X) = ',num2str(VX),])
disp(['Fraction of samples with one or more matches   pm = ', num2str(pm),])
disp(['Probability of one or more matches in a sample Pm = ', num2str(1-P0),])

Distribuciones

comb.m

Descripción del Código:

función comb.m y = comb (n, k) Calcula los coeficientes binomiales. $k$ puede ser una matriz de números enteros entre 0 y $n$ . El resultado $y$ es una matriz de las mismas dimensiones.

Contestar

function y = comb(n,k)
% COMB y = comb(n,k) Binomial coefficients
% Version of 12/10/92
% Computes binomial coefficients C(n,k)
% k may be a matrix of integers between 0 and n
% result y is a matrix of the same dimensions
y = round(gamma(n+1)./(gamma(k + 1).*gamma(n + 1 - k)));

ibinom.m

Descripción del Código:

ibinom.m Distribución binomial — términos individuales. Tenemos dos funciones m ibinom y cbinom para calcular términos individuales y acumulativos $P(S_n = k)$ y $P(S_n \ge k)$ , respectivamente.

$P(S_n = k) = C(n, k) p^k (1 - p)^{n - k}$ y $P(S_n \ge k) = \sum_{r = k}^{n} P(S_n = r)$ $0 \le k \le n$

Para estas funciones m, se utiliza una modificación de una estrategia de cómputo empleada por S. Weintraub: Tablas de la Distribución de Probabilidad Binomial Acumulada para Valores Pequeños de p, 1963. El libro contiene un análisis de errores particularmente útil, escrito por Leo J. Cohen. La experimentación con sumas y expectativas indica una precisión para los cálculos de ibinom y cbinom que es mejor que $10^{-10}$ para $n = 1000$ y $p$ de 0.01 a 0.99. Una precisión similar se mantiene para valores de $n$ hasta 5000, siempre $np$ o $nq$ están limitados a aproximadamente 500. Por encima de este valor para $np$ o $nq$ , los cálculos se descomponen. Para términos individuales, la función y = ibinom (n, p, k) calcula las probabilidades para $n$ un entero positivo, $k$ una matriz de enteros entre 0 y $n$ . La salida es una matriz de las probabilidades binomiales correspondientes.

Contestar

function y = ibinom(n,p,k)
% IBINOM  y = ibinom(n,p,k) Individual binomial probabilities
% Version of 10/5/93
% n is a positive integer; p is a probability
% k a matrix of integers between 0 and n
% y = P(X>=k) (a matrix of probabilities)
if p > 0.5
a = [1 ((1-p)/p)*ones(1,n)];
b = [1 n:-1:1];
c = [1 1:n];
br = (p^n)*cumprod(a.*b./c);
bi = fliplr(br);
else
a = [1 (p/(1-p))*ones(1,n)];
b = [1 n:-1:1];
c = [1 1:n];
bi = ((1-p)^n)*cumprod(a.*b./c);
end
y = bi(k+1);

ipoisson.m

Descripción del Código:

ipoisson.m Distribución de Poisson — términos individuales. Como en el caso de la distribución binomial, tenemos una función m para los términos individuales y otra para el caso acumulativo. Las funciones m ipoisson y cpoisson utilizan una estrategia computacional similar a la utilizada para el caso binomial. Esto no sólo funciona para grandes $\mu$ , sino que la precisión es al menos tan buena como la de las funciones binomiales m. La experiencia indica que las funciones m son buenas para $\mu \le 700$ . Se descomponen en aproximadamente 710, en gran parte debido a las limitaciones de la función exponencial de MATLAB. Para términos individuales, la función y = ipoisson (mu, k) calcula las probabilidades para $mu$ un número entero positivo, $k$ un vector de fila o columna de enteros no negativos. La salida es un vector de fila de las probabilidades de Poisson correspondientes.

Contestar

function y = ipoisson(mu,k)
% IPOISSON y = ipoisson(mu,k) Individual Poisson probabilities
% Version of 10/15/93
% mu = mean value
% k may be a row or column vector of integer values
% y = P(X = k) (a row vector of probabilities)
K = max(k);
p = exp(-mu)*cumprod([1 mu*ones(1,K)]./[1 1:K]);
y = p(k+1);

cpoisson.m

Descripción del Código:

cpoisson.m Distribución de Poisson: términos acumulativos. función y = cpoisson (mu, k), calcula $P(X \ge k)$ , donde $k$ puede haber una fila o un vector de columna de enteros no negativos. La salida es un vector de fila de las probabilidades correspondientes.

Contestar

function y = cpoisson(mu,k)
% CPOISSON y = cpoisson(mu,k) Cumulative Poisson probabilities
% Version of 10/15/93
% mu = mean value mu 
% k may be a row or column vector of integer values
% y = P(X >= k) (a row vector of probabilities)
K = max(k);
p = exp(-mu)*cumprod([1 mu*ones(1,K)]./[1 1:K]);
pc = [1 1 - cumsum(p)];
y = pc(k+1);

nbinom.m

Descripción del Código:

nbinom.m Binomial negativo — función y = nbinom (m, p, k) calcula la probabilidad de que el éxito $m$ th en una secuencia de Bernoulli ocurra en el $k$ th ensayo.

Contestar

function y = nbinom(m, p, k)
% NBINOM y = nbinom(m, p, k) Negative binomial probabilities
% Version of 12/10/92
% Probability the mth success occurs on the kth trial
% m a positive integer;  p a probability
% k a matrix of integers greater than or equal to m
% y = P(X=k) (a matrix of the same dimensions as k)
q = 1 - p;
y = ((p^m)/gamma(m)).*(q.^(k - m)).*gamma(k)./gamma(k - m + 1);

gaussian.m

Descitación de Código:

gaussian.m función y = gaussiana (m, v, t) calcula la función de distribución gaussiana (Normal) para el valor medio $m$ $v$ , varianza y matriz $t$ de valores. El resultado $y = P(X \le t)$ es una matriz de las mismas dimensiones que $t$ .

Contestar

function y = gaussian(m,v,t)
% GAUSSIAN y = gaussian(m,v,t) Gaussian distribution function
% Version of 11/18/92
% Distribution function for X ~ N(m, v)
% m = mean,  v = variance
% t is a matrix of evaluation points
% y = P(X<=t) (a matrix of the same dimensions as t)
u = (t - m)./sqrt(2*v);
if u >= 0
        y = 0.5*(erf(u) + 1);
else
        y = 0.5*erfc(-u);
end

gaussdensity.m

Descripción del Código:

gaussdensity.m función y = gaussdensidad (m, v, t) calcula la función de densidad gaussiana $f_X (t)$ para el valor medio $m$ $t$ , varianza y matriz $t$ de valores.

Contestar

function y = gaussdensity(m,v,t)
% GAUSSDENSITY y = gaussdensity(m,v,t) Gaussian density
% Version of 2/8/96
% m = mean,  v = variance
% t is a matrix of evaluation points
y = exp(-((t-m).^2)/(2*v))/sqrt(v*2*pi);

norminv.m

Descripción del Código:

función norminv.m y = norminv (m, v, p) calcula la inversa (la función cuantil) de la función de distribución gaussiana para el valor medio $m$ $v$ , varianza y $p$ una matriz de probabilidades.

Contestar

function y = norminv(m,v,p)
% NORMINV y = norminv(m,v,p) Inverse gaussian distribution
% (quantile function for gaussian)
% Version of 8/17/94
% m = mean,  v = variance
% t is a matrix of evaluation points
if p >= 0
  u = sqrt(2)*erfinv(2*p - 1);
else
  u = -sqrt(2)*erfinv(1 - 2*p);
end
y = sqrt(v)*u + m;

gammadbn.m

Descripción del Código:

función gammadbn.m y = gammadbn (alpha, lambda, t) calcula la función de distribución para una distribución gamma con parámetros alpha, lambda. $t$ es una matriz de puntos de evaluación. El resultado es una matriz del mismo tamaño.

Contestar

function y = gammadbn(alpha, lambda, t)
% GAMMADBN y = gammadbn(alpha, lambda, t) Gamma distribution
% Version of 12/10/92
% Distribution function for X ~ gamma (alpha, lambda)
% alpha, lambda are positive parameters
% t may be a matrix of positive numbers
% y = P(X<= t) (a matrix of the same dimensions as t)
y = gammainc(lambda*t, alpha);

beta.m

Descripción del Código:

función beta.m y = beta (r, s, t) calcula la función de densidad para la distribución beta con parámetros $r, s, t$ es una matriz de números entre cero y uno. El resultado es una matriz del mismo tamaño.

Contestar

function y = beta(r,s,t)
% BETA y = beta(r,s,t) Beta density function
% Version of 8/5/93
% Density function for Beta (r,s) distribution
% t is a matrix of evaluation points between 0 and 1
% y is a matrix of the same dimensions as t
y = (gamma(r+s)/(gamma(r)*gamma(s)))*(t.^(r-1).*(1-t).^(s-1));

betadbn.m

Descripción del Código:

betadbn.m función y = betadbn (r, s, t) calcula la función de distribución para la distribución beta con parámetros $r, s, t$ es una matriz de puntos de evaluación. El resultado es una matriz del mismo tamaño.

Contestar

function y = betadbn(r,s,t)
% BETADBN y = betadbn(r,s,t) Beta distribution function
% Version of 7/27/93
% Distribution function for X  beta(r,s)
% y = P(X<=t) (a matrix of the same dimensions as t)
y = betainc(t,r,s);

weibull.m

Descripción del Código:

weibull.m función y = weibull (alpha, lambda, t) calcula la función de densidad para la distribución Weibull con parámetros alpha, lambda. $t$ es una matriz de puntos de evaluación. El resultado es una matriz del mismo tamaño.

Contestar

function y = weibull(alpha,lambda,t)
% WEIBULL y = weibull(alpha,lambda,t) Weibull density 
% Version of 1/24/91
% Density function for X ~ Weibull (alpha, lambda, 0)
% t is a matrix of positive evaluation points
% y is a matrix of the same dimensions as t
y = alpha*lambda*(t.^(alpha - 1)).*exp(-lambda*(t.^alpha));

weibulld.m

Descripción del Código:

weibulld.m función y = weibulld (alpha, lambda, t) calcula la función de distribución para la distribución Weibull con parámetros alpha, lambda. $t$ es una matriz de puntos de evaluación. El resultado es una matriz del mismo tamaño.

Contestar

function y = weibulld(alpha, lambda, t)
% WEIBULLD y = weibulld(alpha, lambda, t) Weibull distribution function
% Version of 1/24/91
% Distribution function for X ~ Weibull (alpha, lambda, 0)
% t is a matrix of positive evaluation points
% y = P(X<=t) (a matrix of the same dimensions as t)
y = 1 - exp(-lambda*(t.^alpha));

Dstribuciones binomiales, de Poisson y Gaussianas

bincomp.m

Descripción del Código:

bincomp.m Comparación gráfica de las distribuciones binomiales, Poisson y Gaussianas. El procedimiento requiere parámetros binomiales $n, p$ , determina un rango razonable de puntos de evaluación y traza en la misma gráfica la función de distribución binomial, la función de distribución de Poisson y la función de distribución gaussiana con el ajuste llamado “corrección de continuidad”.

Contestar

% BINCOMP file bincomp.m  Approx of binomial by Poisson and gaussian
% Version of 5/24/96
% Gaussian adjusted for "continuity correction"
% Plots distribution functions for specified parameters n, p
n = input('Enter the parameter n  ');
p = input('Enter the parameter p  ');
a = floor(n*p-2*sqrt(n*p));
a = max(a,1);                         % Prevents zero or negative indices
b = floor(n*p+2*sqrt(n*p));
k = a:b; 
Fb = cumsum(ibinom(n,p,0:n));         % Binomial distribution function
Fp = cumsum(ipoisson(n*p,0:n));       % Poisson distribution function
Fg = gaussian(n*p,n*p*(1 - p),k+0.5); % Gaussian distribution function
stairs(k,Fb(k+1))                     % Plotting details
hold on
plot(k,Fp(k+1),'-.',k,Fg,'o') 
hold off
xlabel('t values')                    % Graph labeling details
ylabel('Distribution function')
title('Approximation of Binomial by Poisson and Gaussian')
grid 
legend('Binomial','Poisson','Adjusted Gaussian')
disp('See Figure for results')

poissapp.m

Descitación de Código:

poissapp.m Comparación gráfica de las distribuciones de Poisson y Gauss. El procedimiento requiere un valor del parámetro de Poisson mu, luego calcula y traza la función de distribución de Poisson, la función de distribución gaussiana y la función de distribución gaussiana ajustada.

Contestar

% POISSAPP file poissapp.m  Comparison of Poisson and gaussian
% Version of 5/24/96
% Plots distribution functions for specified parameter mu
mu = input('Enter the parameter mu  ');
n = floor(1.5*mu);
k = floor(mu-2*sqrt(mu)):floor(mu+2*sqrt(mu));              
FP = cumsum(ipoisson(mu,0:n));
FG = gaussian(mu,mu,k); 
FC = gaussian(mu,mu,k-0.5);      
stairs(k,FP(k))                 
hold on
plot(k,FG,'-.',k,FC,'o') 
hold off
grid
xlabel('t values')
ylabel('Distribution function')
title('Gaussian Approximation to Poisson Distribution')
legend('Poisson','Gaussian','Adjusted Gaussian')
disp('See Figure for results')

Configuración para variables aleatorias simples

Si una variable aleatoria simple $X$ está en forma canónica, la distribución consiste en los coeficientes de las funciones indicadoras (los valores de $X$ ) y las probabilidades de los eventos correspondientes. Si $X$ está en una forma primitiva distinta a la canónica, la operación csort se aplica a los coeficientes de las funciones indicadoras y las probabilidades de los eventos correspondientes para obtener la distribución. Si $Z = g(X)$ y $X$ está en una forma primitiva, entonces el valor de $Z$ on el evento en la partición asociada con $t_i$ es $g(t_i)$ . La distribución para Z se obtiene aplicando csort a la $g(t_i)$ y la $p_i$ . De igual manera, si $Z = g(X, Y)$ y la distribución conjunta está disponible, el valor $g(t_i, u_j)$ está asociado con $P(X = t_i, Y = u_j)$ . La distribución para $Z$ se obtiene aplicando csort a la matriz de valores y a la matriz de probabilidades correspondiente.

canónica.m

Descripción del Código:

canonic.m El procedimiento determina la distribución para una variable aleatoria simple en forma afín, cuando las probabilidades minterm están disponibles. Los datos de entrada son un vector de fila de coeficientes para las funciones indicadoras en la forma afín (con el valor constante último) y un vector de fila de las probabilidades del minterm generado por los eventos. Los resultados consisten en un vector de fila de valores y un vector de fila de las probabilidades correspondientes.

Contestar

% CANONIC file canonic.m Distribution for simple rv in affine form
% Version of 6/12/95
% Determines the distribution for a simple random variable
% in affine form, when the minterm probabilities are available.
% Uses the m-functions mintable and csort.
% The coefficient vector must contain the constant term.
% If the constant term is zero, enter 0 in the last place.
c  = input(' Enter row vector of coefficients  ');
pm = input(' Enter row vector of minterm probabilities  ');
n  = length(c) - 1;
if 2^n ~= length(pm)
   error('Incorrect minterm probability vector length');
end
M  = mintable(n);            % Provides a table of minterm patterns
s  = c(1:n)*M + c(n+1);      % Evaluates X on each minterm
[X,PX] = csort(s,pm);        % s = values; pm = minterm probabilities
XDBN = [X;PX]';
disp('Use row matrices X and PX for calculations')
disp('Call for XDBN to view the distribution')

canonicf.m

Descripción del Código:

canonicf.m function [x, px] = canonicf (c, pm) es una versión de función de canonic, que permite nombrar arbitrariamente las variables.

Contestar

function [x,px] = canonicf(c,pm)
% CANONICF  [x,px] = canonicf(c,pm)  Function version of canonic
% Version of 6/12/95
% Allows arbitrary naming of variables
n = length(c) - 1;
if 2^n ~= length(pm)
   error('Incorrect minterm probability vector length');
end
M  = mintable(n);              % Provides a table of minterm patterns
s  = c(1:n)*M + c(n+1);        % Evaluates X on each minterm
[x,px]  = csort(s,pm);         % s = values; pm = minterm probabilities

jcalc.m

Descripción del Código:

jcalc.m Se configura para cálculos para variables aleatorias simples conjuntas. La matriz $P$ de $P(X = t_i, Y = u_j)$ está dispuesta como en el plano (es decir, valores de $Y$ aumento hacia arriba). La función de MATLAB meshgrid se aplica a la matriz de fila $X$ y $Y$ a la matriz de fila invertida para poner un $X$ valor y un $Y$ valor apropiado en cada posición. Estas se encuentran en las “matrices de cálculo” $t$ y $u$ , respectivamente, las cuales se utilizan en la determinación de probabilidades y expectativas de diversas funciones de $t, u$ .

Contestar

% JCALC  file jcalc.m  Calculation setup for joint simple rv
% Version of 4/7/95 (Update of prompt and display 5/1/95)
% Setup for calculations for joint simple random variables
% The joint probabilities arranged as on the plane 
% (top row corresponds to largest value of Y) 
P = input('Enter JOINT PROBABILITIES (as on the plane)  ');
X = input('Enter row matrix of VALUES of X  ');
Y = input('Enter row matrix of VALUES of Y  ');
PX = sum(P);            % probabilities for X
PY = fliplr(sum(P'));   % probabilities for Y
[t,u] = meshgrid(X,fliplr(Y));
disp(' Use array operations on matrices X, Y, PX, PY, t, u, and P')

jcalcf.m

Descripción del Código:

jcalcf.m función [x, y, t, u, px, py, p] = jcalcf (X, Y, P) es una versión de función de jcalc, que permite nombrar arbitrariamente las variables.

Contestar

function [x,y,t,u,px,py,p] = jcalcf(X,Y,P)
% JCALCF [x,y,t,u,px,py,p] = jcalcf(X,Y,P) Function version of jcalc
% Version of 5/3/95
% Allows arbitrary naming of variables
if sum(size(P) ~= [length(Y) length(X)]) > 0
  error('     Incompatible vector sizes')
end
x = X;
y = Y;
p = P;
px = sum(P);
py = fliplr(sum(P'));
[t,u] = meshgrid(X,fliplr(Y));

jointzw.m

Descripción del Código:

jointzw.m Establece la distribución conjunta para $Z = g(X, Y)$ y $W = h(X, Y)$ y proporciona matrices de cálculo como en jcalc. Las entradas son $P, X$ y así $Y$ como expresiones de matriz para $g(t, u)$ y $h(t, u)$ . Las salidas son matrices $Z, W, PZW$ para la distribución conjunta $PZ, PW$ , las probabilidades marginales y las matrices de cálculo $v, w$ .

Contestar

% JOINTZW file jointzw.m Joint dbn for two functions of (X,Y)
% Version of 4/29/97
% Obtains joint distribution for
% Z = g(X,Y) and W = h(X,Y)
% Inputs P, X, and Y as well as array
% expressions for g(t,u) and h(t,u)
P = input('Enter joint prob for (X,Y) ');
X = input('Enter values for X ');
Y = input('Enter values for Y ');
[t,u] = meshgrid(X,fliplr(Y));
G = input('Enter expression for g(t,u) ');
H = input('Enter expression for h(t,u) ');
[Z,PZ] = csort(G,P);
[W,PW] = csort(H,P);
r = length(W);
c = length(Z);
PZW = zeros(r,c);
for i = 1:r
  for j = 1:c
   a = find((G==Z(j))&(H==W(i)));
   if ~isempty(a)
    PZW(i,j) = total(P(a));
   end
  end
end
PZW = flipud(PZW);
[v,w] = meshgrid(Z,fliplr(W));
if (G==t)&(H==u)
  disp(' ')
  disp('  Note:  Z = X and W = Y')
  disp(' ')
elseif  G==t
  disp(' ')
  disp('  Note:  Z = X')
  disp(' ')
elseif H==u
  disp(' ')
  disp('  Note:  W = Y')
  disp(' ')
end
disp('Use array operations on Z, W, PZ, PW, v, w, PZW')

jdtest.m

Descripción del Código:

jdtest.m Prueba una matriz de probabilidad conjunta $P$ para entradas negativas y probabilidad total unitaria..

Contestar

function y = jdtest(P)
% JDTEST y = jdtest(P) Tests P for unit total and negative elements
% Version of 10/8/93
M = min(min(P));
S = sum(sum(P));
if M < 0
  y = 'Negative entries';
elseif abs(1 - S) > 1e-7
  y = 'Probabilities do not sum to one';
else
  y = 'P is a valid distribution';
end

Configuración para variables aleatorias generales

tappr.m

Descripción del Código:

tappr.m Utiliza la función de densidad para establecer una aproximación discreta a la distribución para una variable aleatoria absolutamente continua $X$ .

Contestar

% TAPPR file tappr.m  Discrete approximation to ac random variable
% Version of 4/16/94
% Sets up discrete approximation to distribution for
% absolutely continuous random variable  X
% Density is entered as a function of t
r = input('Enter matrix [a b] of x-range endpoints  ');
n = input('Enter number of x approximation points  ');
d = (r(2) - r(1))/n;
t = (r(1):d:r(2)-d) +d/2;
PX = input('Enter density as a function of t  ');
PX = PX*d;
PX = PX/sum(PX);
X  = t;
disp('Use row matrices X and PX as in the simple case')

tuappr.m

Descripción del Código:

tuappr.m Utiliza la densidad de juntas para establecer aproximaciones discretas a $X, Y, t, u$ , y densidad.

Contestar

% TUAPPR file tuappr.m Discrete approximation to joint ac pair
% Version of 2/20/96
% Joint density entered as a function of t, u
% Sets up discrete approximations to X, Y, t, u, and density
rx = input('Enter matrix [a b] of X-range endpoints  ');
ry = input('Enter matrix [c d] of Y-range endpoints  ');
nx = input('Enter number of X approximation points  ');
ny = input('Enter number of Y approximation points  ');
dx = (rx(2) - rx(1))/nx;
dy = (ry(2) - ry(1))/ny;
X  = (rx(1):dx:rx(2)-dx) + dx/2;
Y  = (ry(1):dy:ry(2)-dy) + dy/2;
[t,u] = meshgrid(X,fliplr(Y));
P  = input('Enter expression for joint density  ');
P  = dx*dy*P;
P  = P/sum(sum(P));
PX = sum(P);
PY = fliplr(sum(P'));
disp('Use array operations on X, Y, PX, PY, t, u, and P')

dfappr.m

dfappr.m Distribución discreta aproximada de la función de distribución ingresada en función de $t$ .

Contestar

% DFAPPR file dfappr.m Discrete approximation from distribution function
% Version of 10/21/95
% Approximate discrete distribution from distribution
% function entered as a function of t
r = input('Enter matrix [a b] of X-range endpoints  ');
s = input('Enter number of X approximation points  ');
d = (r(2) - r(1))/s;
t = (r(1):d:r(2)-d) +d/2;
m  = length(t);
f  = input('Enter distribution function F as function of t  ');
f  = [0 f];
PX = f(2:m+1) - f(1:m);
PX = PX/sum(PX);
X  = t - d/2;
disp('Distribution is in row matrices X and PX')

acsetup.m

Descripción del Código:

acsetup.m Distribución aproximada para variable aleatoria absolutamente continua $X$ . Densidad se ingresa como una función variable de cadena de $t$ .

Contestar

% ACSETUP file acsetup.m Discrete approx from density as string variable
% Version of 10/22/94
% Approximate distribution for absolutely continuous rv X
% Density is entered as a string variable function of t
disp('DENSITY f is entered as a STRING VARIABLE.')
disp('either defined previously or upon call.')
r  = input('Enter matrix [a b] of x-range endpoints  ');
s  = input('Enter number of x approximation points  ');
d  = (r(2) - r(1))/s;
t  = (r(1):d:r(2)-d) +d/2;
m  = length(t);
f  = input('Enter density as a function of t  ');
PX = eval(f);
PX = PX*d;
PX = PX/sum(PX);
X  = t;
disp('Distribution is in row matrices X and PX')

dfsetup.m

Descripción del Código:

dfsetup.m Distribución discreta aproximada de la función de distribución ingresada como una función variable de cadena de $t$ .

Contestar

% DFSETUP file dfsetup.m  Discrete approx from string dbn function
% Version of 10/21/95
% Approximate discrete distribution from distribution
% function entered as string variable function of t
disp('DISTRIBUTION FUNCTION F is entered as a STRING')
disp('VARIABLE, either defined previously or upon call')
r = input('Enter matrix [a b] of X-range endpoints  ');
s = input('Enter number of X approximation points  ');
d = (r(2) - r(1))/s;
t = (r(1):d:r(2)-d) +d/2;
m  = length(t);
F  = input('Enter distribution function F as function of t  ');
f  = eval(F);
f  = [0 f];
PX = f(2:m+1) - f(1:m);
PX = PX/sum(PX);
X  = t - d/2;
disp('Distribution is in row matrices X and PX')

Configuración para variables aleatorias simples independientes

MATLAB versión 5.1 cuenta con provisiones para arreglos multidimensionales, que hacen posible una implementación más directa de icalc3 e icalc4.

icalc.m

Descripción del Código:

icalc.m Configuración de cálculo para un par independiente de variables aleatorias simples. La entrada consiste en distribuciones marginales para $X, Y$ . La salida es la distribución conjunta y el cálculo de matrices $t, u$ .

Contestar

% ICALC file icalc.m  Calculation setup for independent pair
% Version of 5/3/95
% Joint calculation setup for independent pair
X  = input('Enter row matrix of X-values  ');
Y  = input('Enter row matrix of Y-values  ');
PX = input('Enter X probabilities  ');
PY = input('Enter Y probabilities  ');
[a,b] = meshgrid(PX,fliplr(PY));
P  = a.*b;                      % Matrix of joint independent probabilities 
[t,u] = meshgrid(X,fliplr(Y));  % t, u matrices for joint calculations
disp(' Use array operations on matrices X, Y, PX, PY, t, u, and P')

icalcf.m

icalcf.m [x, y, t, u, px, py, p] = icalcf (X, Y, PX, PY) es una versión de función de icalc, que permite nombrar arbitrariamente las variables.

Contestar

function [x,y,t,u,px,py,p] = icalcf(X,Y,PX,PY)
% ICALCF [x,y,t,u,px,py,p] = icalcf(X,Y,PX,PY) Function version of icalc
% Version of 5/3/95
% Allows arbitrary naming of variables
x = X;
y = Y;
px = PX;
py = PY;
if length(X) ~= length(PX)
  error('     X and PX of different lengths')
elseif length(Y) ~= length(PY)
  error('     Y and PY of different lengths')
end
[a,b] = meshgrid(PX,fliplr(PY));
p   = a.*b;                       % Matrix of joint independent probabilities 
[t,u] = meshgrid(X,fliplr(Y));    % t, u matrices for joint calculations

icalc3.m

Descripción del Código:

icalc3.m Configuración de cálculo para una clase independiente de tres variables aleatorias simples.

Contestar

% ICALC3 file icalc3.m Setup for three independent rv
% Version of 5/15/96
% Sets up for calculations for three
% independent simple random variables
% Uses m-functions rep, elrep, kronf
X  = input('Enter row matrix of X-values  ');
Y  = input('Enter row matrix of Y-values  ');
Z  = input('Enter row matrix of Z-values  ');
PX = input('Enter X probabilities  ');
PY = input('Enter Y probabilities  ');
PZ = input('Enter Z probabilities  ');
n  = length(X);
m  = length(Y);
s  = length(Z);
[t,u] = meshgrid(X,Y);
t  = rep(t,1,s);
u  = rep(u,1,s);
v  = elrep(Z,m,n);  % t,u,v matrices for joint calculations
P  = kronf(PZ,kronf(PX,PY'));
disp('Use array operations on matrices X, Y, Z,')
disp('PX, PY, PZ, t, u, v, and P')

icalc4.m

Descripción del Código:

icalc4.m Configuración de cálculo para una clase independiente de cuatro variables aleatorias simples.

Contestar

% ICALC4 file icalc4.m Setup for four independent rv
% Version of 5/15/96
% Sets up for calculations for four
% independent simple random variables
% Uses m-functions rep, elrep, kronf
X  = input('Enter row matrix of X-values  ');
Y  = input('Enter row matrix of Y-values  ');
Z  = input('Enter row matrix of Z-values  ');
W  = input('Enter row matrix of W-values  ');
PX = input('Enter X probabilities  ');
PY = input('Enter Y probabilities  ');
PZ = input('Enter Z probabilities  ');
PW = input('Enter W probabilities  ');
n  = length(X);
m  = length(Y);
s  = length(Z);
r  = length(W);
[t,u] = meshgrid(X,Y);
t = rep(t,r,s);
u = rep(u,r,s);
[v,w] = meshgrid(Z,W);
v = elrep(v,m,n);  % t,u,v,w matrices for joint calculations
w = elrep(w,m,n);
P = kronf(kronf(PZ,PW'),kronf(PX,PY'));
disp('Use array operations on matrices X, Y, Z, W')
disp('PX, PY, PZ, PW, t, u, v, w, and P')

Cálculos para variables aleatorias

ddbn.m

Descripción del Código:

ddbn.m Utiliza la distribución de una variable aleatoria simple (o aproximación simple) para trazar un gráfico de pasos para la función de distribución $F_X$

Contestar

% DDBN file ddbn.m Step graph of distribution function
% Version of 10/25/95
% Plots step graph of dbn function FX from 
% distribution of simple rv (or simple approximation)
xc  = input('Enter row matrix of VALUES  ');
pc = input('Enter row matrix of PROBABILITIES  ');
m  = length(xc);
FX = cumsum(pc);
xt = [xc(1)-1-0.1*abs(xc(1)) xc xc(m)+1+0.1*abs(xc(m))];
FX = [0 FX 1];        % Artificial extension of range and domain
stairs(xt,FX)         % Plot of stairstep graph
hold on
plot(xt,FX,'o')       % Marks values at jump
hold off
grid                  
xlabel('t')
ylabel('u = F(t)')
title('Distribution Function')

cdbn.m

Descripción del Código:

cdbn.m Traza una gráfica continua de una función de distribución de una variable aleatoria simple (o aproximación simple).

Contestar

% CDBN file cdbn.m Continuous graph of distribution function
% Version of 1/29/97
% Plots continuous graph of dbn function FX from 
% distribution of simple rv (or simple approximation)
xc  = input('Enter row matrix of VALUES  ');
pc = input('Enter row matrix of PROBABILITIES  ');
m  = length(xc);
FX = cumsum(pc);
xt = [xc(1)-0.01 xc xc(m)+0.01];
FX = [0 FX FX(m)];    % Artificial extension of range and domain
plot(xt,FX)           % Plot of continuous graph
grid                  
xlabel('t')
ylabel('u = F(t)')
title('Distribution Function')

simple.m

Descripción del Código:

simple.m Calcula cuantitos básicos para variables aleatorias simples a partir de la distribución, entrada como matrices de fila $X$ y $PX$ .

Contestar

% SIMPLE file simple.m Calculates basic quantites for simple rv
% Version of 6/18/95
X  = input('Enter row matrix of X-values  ');
PX = input('Enter row matrix PX of X probabilities  ');
n  = length(X);          % dimension of X
EX = dot(X,PX)           % E[X]
EX2 = dot(X.^2,PX)       % E[X^2]
VX = EX2 - EX^2          % Var[X]
disp(' ')
disp('Use row matrices X and PX for further calculations')

jddbn.m

Descripción del Código:

jddbn.m Representación de la función de distribución conjunta para par simple mediante la obtención del valor de $F_{XY}$ en las esquinas inferiores izquierdas de cada celda de cuadrícula.

Contestar

% JDDBN file jddbn.m  Joint distribution function
% Version of 10/7/96
% Joint discrete distribution function for
% joint  matrix P (arranged as on the plane).
% Values at lower left hand corners of grid cells
P = input('Enter joint probability matrix (as on the plane)  ');
FXY = flipud(cumsum(flipud(P)));
FXY = cumsum(FXY')';
disp('To view corner values for joint dbn function, call for FXY')

jsimple.m

Descripción del Código:

jsimple.m Calcula cantidades básicas para un par simple conjunto $\{X, Y\}$ a partir de la distribución articular $X, Y, P$ como en jcalc. Las cantidades calculadas incluyen medias, varianzas, covarianza, línea de regresión y curva de regresión (expectativa condicional $E[Y|X = t]$

Contestar

% JSIMPLE file jsimple.m  Calculates basic quantities for joint simple rv
% Version of 5/25/95
% The joint probabilities are arranged as on the plane 
% (the top row corresponds to the largest value of Y) 
P = input('Enter JOINT PROBABILITIES (as on the plane)  ');
X = input('Enter row matrix of VALUES of X  ');
Y = input('Enter row matrix of VALUES of Y  ');
disp(' ')
PX = sum(P);               % marginal distribution for X
PY = fliplr(sum(P'));      % marginal distribution for Y
XDBN = [X; PX]';
YDBN = [Y; PY]';
PT  = idbn(PX,PY);
D  = total(abs(P - PT));   % test for difference
if D > 1e-8                % to prevent roundoff error masking zero
  disp('{X,Y} is NOT independent')
 else
  disp('{X,Y} is independent')
end
disp(' ')
[t,u] = meshgrid(X,fliplr(Y));
EX  = total(t.*P)          % E[X]
EY  = total(u.*P)          % E[Y]
EX2 = total((t.^2).*P)     % E[X^2]
EY2 = total((u.^2).*P)     % E[Y^2]
EXY = total(t.*u.*P)       % E[XY]
VX  = EX2 - EX^2           % Var[X]
VY  = EY2 - EY^2           % Var[Y]
cv = EXY - EX*EY;          % Cov[X,Y] = E[XY] - E[X]E[Y]
if abs(cv) > 1e-9          % to prevent roundoff error masking zero
   CV = cv
 else
   CV = 0
end
a = CV/VX                  % regression line of Y on X is
b = EY - a*EX              %       u = at + b
R = CV/sqrt(VX*VY);        % correlation coefficient rho
disp(['The regression line of Y on X is:  u = ',num2str(a),'t + ',num2str(b),])
disp(['The correlation coefficient is:  rho = ',num2str(R),])
disp(' ')
eYx = sum(u.*P)./PX; 
EYX = [X;eYx]';
disp('Marginal dbns are in X, PX, Y, PY; to view, call XDBN, YDBN')
disp('E[Y|X = x] is in eYx; to view, call for EYX')
disp('Use array operations on matrices X, Y, PX, PY, t, u, and P')

japprox.m

Descripción del Código:

japprox.m Asume configuración discreta y calcula cantidades básicas para un par de variables aleatorias como en jsimple. Traza la línea de regresión y la curva de regresión.

Contestar

% JAPPROX file japprox.m Basic quantities for ac pair {X,Y}
% Version of 5/7/96
% Assumes tuappr has set X, Y, PX, PY, t, u, P
EX  = total(t.*P)          % E[X]
EY  = total(u.*P)          % E[Y]
EX2 = total(t.^2.*P)       % E[X^2]
EY2 = total(u.^2.*P)       % E[Y^2]
EXY = total(t.*u.*P)       % E[XY]
VX  = EX2 - EX^2           % Var[X]
VY  = EY2 - EY^2           % Var[Y]
cv = EXY - EX*EY;          % Cov[X,Y] = E[XY] - E[X]E[Y]
if abs(cv) > 1e-9  % to prevent roundoff error masking zero
   CV = cv
 else
   CV = 0
end
a = CV/VX                  % regression line of Y on X is
b = EY - a*EX              % u = at + b
R = CV/sqrt(VX*VY);
disp(['The regression line of Y on X is:  u = ',num2str(a),'t + ',num2str(b),])
disp(['The correlation coefficient is:  rho = ',num2str(R),])
disp(' ')
eY = sum(u.*P)./sum(P);    % eY(t) = E[Y|X = t]
RL = a*X + b;
plot(X,RL,X,eY,'-.')
grid
title('Regression line and Regression curve')
xlabel('X values')
ylabel('Y values')
legend('Regression line','Regression curve')
clear eY                   % To conserve memory
clear RL
disp('Calculate with X, Y, t, u, P, as in joint simple case')

Cálculos y pruebas para variables aleatorias independientes

mgsum.m

Descripción del Código:

mgsum.m function [z, pz] = mgsum (x, y, px, py) determina la distribución para la suma de un par independiente de variables aleatorias simples a partir de sus distribuciones.

Contestar

function [z,pz] = mgsum(x,y,px,py)
% MGSUM [z,pz] = mgsum(x,y,px,py)  Sum of two independent simple rv
% Version of 5/6/96
% Distribution for the sum of two independent simple random variables
% x is a vector (row or column) of X values  
% y is a vector (row or column) of Y values
% px is a vector (row or column) of X probabilities
% py is a vector (row or column) of Y probabilities
% z and pz are row vectors
[a,b] = meshgrid(x,y);
t  = a+b;
[c,d] = meshgrid(px,py);
p  = c.*d;
[z,pz]  = csort(t,p);

mgsum3.m

Descripción del Código:

mgsum3.m función [w, pw] = mgsum3 (x, y, z, px, py, pz) extiende mgsum a tres variables aleatorias mediante la aplicación repetida de mgsum. De manera similar para mgsum4.m.

Contestar

function [w,pw] = mgsum3(x,y,z,px,py,pz)
% MGSUM3 [w,pw] = mgsum3(x,y,z,px,py,y) Sum of three independent simple rv
% Version of 5/2/96
% Distribution for the sum of three independent simple random variables
% x is a vector (row or column) of X values  
% y is a vector (row or column) of Y values
% z is a vector (row or column) of Z values
% px is a vector (row or column) of X probabilities
% py is a vector (row or column) of Y probabilities
% pz is a vector (row or column) of Z probabilities
% W and pW are row vectors
[a,pa] = mgsum(x,y,px,py);
[w,pw] = mgsum(a,z,pa,pz);

mgnsum.m

Descripción del Código:

función mgnsum.m [z, pz] = mgnsum (X, P) determina la distribución para una suma de variables aleatorias $n$ independientes. $X$ una matriz $n$ -fila de $X$ -valores y $n$ -fila matriz de $P$ -valores (acolchada con ceros, si es necesario, para hacer que todas las filas tengan la misma longitud.

Contestar

function [z,pz] = mgnsum(X,P)
% MGNSUM [z,pz] = mgnsum(X,P)  Sum of n independent simple rv
% Version of 5/16/96
% Distribution for the sum of n independent simple random variables
% X an n-row matrix of X-values
% P an n-row matrix of P-values
% padded with zeros, if necessary
% to make all rows the same length
[n,r] = size(P);
z  = 0;
pz = 1;
for i = 1:n
  x = X(i,:);
  p = P(i,:);
  x = x(find(p>0));
  p = p(find(p>0));
  [z,pz] = mgsum(z,x,pz,p);
end

mgsumn.m

Descripción del Código:

mgsumn.m function [z, pz] = mgsumn (varargin) es una alternativa a mgnsum, utilizando varargin en la versión 5.1 de MATLAB. El llamado es de la forma [z, pz] = mgsumn ([x1; p1], [x2; p2],..., [xn; pn]).

Contestar

function [z,pz] = mgsumn(varargin)
% MGSUMN [z,pz] = mgsumn([x1;p1],[x2;p2], ..., [xn;pn])
% Version of 6/2/97 Uses MATLAB version 5.1
% Sum of n independent simple random variables
% Utilizes distributions in the form [x;px] (two rows)
% Iterates mgsum
n  = length(varargin);   % The number of distributions
z  = 0;                  % Initialization
pz = 1;
for i = 1:n              % Repeated use of mgsum
  [z,pz] = mgsum(z,varargin{i}(1,:),pz,varargin{i}(2,:));
end

diidsum.m

Descripción del Código:

función diidsum.m [x, px] = diidsum (X, PX, n) determina la suma de variables aleatorias simples $n$ iid, con la distribución común $X$ , $PX$

Contestar

function [x,px] = diidsum(X,PX,n)
% DIIDSUM [x,px] = diidsum(X,PX,n) Sum of n iid simple random variables
% Version of 10/14/95 Input rev 5/13/97
% Sum of n iid rv with common distribution X, PX
% Uses m-function mgsum
x  = X;                       % Initialization
px = PX;
for i = 1:n-1
  [x,px] = mgsum(x,X,px,PX);
end

itest.m

Descripción del Código:

itest.m Prueba para la independencia la matriz $P$ de probabilidades conjuntas para un par simple $\{X, Y\}$ de variables aleatorias.

Contestar

% ITEST file itest.m  Tests P for independence
% Version of 5/9/95
% Tests for independence the matrix of joint 
% probabilities for a simple pair {X,Y}
pt = input('Enter matrix of joint probabilities  ');
disp(' ')
px = sum(pt);                  % Marginal probabilities for X
py = sum(pt');                 % Marginal probabilities for Y (reversed)
[a,b] = meshgrid(px,py); 
PT = a.*b;                     % Joint independent probabilities
D  = abs(pt - PT) > 1e-9;      % Threshold set above roundoff
if total(D) > 0
  disp('The pair {X,Y} is NOT independent')
  disp('To see where the product rule fails, call for D')
else
  disp('The pair {X,Y} is independent')
end

idbn.m

Descripción del Código:

idbn.m función p = idbn (px, py) utiliza probabilidades marginales para determinar la matriz de probabilidad conjunta (dispuesta como en el plano) para un par independiente de variables aleatorias simples.

Contestar

function p = idbn(px,py)
% IDBN p = idbn(px,py)  Matrix of joint independent probabilities
% Version of 5/9/95
% Determines joint probability matrix for two independent
% simple random variables (arranged as on the plane)
[a,b] = meshgrid(px,fliplr(py));
p = a.*b

isimple.m

Descripción del Código:

isimple.m Toma como entradas las distribuciones marginales para un par independiente $\{X, Y\}$ de variables aleatorias simples. Establece la matriz de probabilidad de distribución conjunta $P$ como en idbn, y forma las matrices de cálculo $t, u$ como en jcalc. Calcula cantidades básicas y pone a disposición matrices $X$ $Y$ , $PX$ , $PY$ , $P$ , $t$ ,, $u$ , para cálculos adicionales.

Contestar

% ISIMPLE file isimple.m  Calculations for independent simple rv
% Version of 5/3/95
X  = input('Enter row matrix of X-values  ');
Y  = input('Enter row matrix of Y-values  ');
PX = input('Enter X probabilities  ');
PY = input('Enter Y probabilities  ');
[a,b] = meshgrid(PX,fliplr(PY));
P  = a.*b;                      % Matrix of joint independent probabilities 
[t,u] = meshgrid(X,fliplr(Y));  % t, u matrices for joint calculations
EX  = dot(X,PX)                 % E[X]
EY  = dot(Y,PY)                 % E[Y]
VX  = dot(X.^2,PX) - EX^2       % Var[X]
VY  = dot(Y.^2,PY) - EY^2       % Var[Y]
disp(' Use array operations on matrices X, Y, PX, PY, t, u, and P')

Funciones cuantiles para distribuciones delimitadas

dquant.m

Descripción del Código:

dquant.m función t = dquant (X, PX, U) determina los valores de la función cuantil para una variable aleatoria simple con distribución $X$ , $PX$ a los valores de probabilidad en el vector de fila $U$ . El vector de probabilidad $U$ suele ser determinado por un generador de números aleatorios.

Contestar

function t = dquant(X,PX,U)
% DQUANT t = dquant(X,PX,U)  Quantile function for a simple random variable
% Version of 10/14/95
% U is a vector of probabilities
m  = length(X);
n  = length(U);
F  = [0 cumsum(PX)+1e-12]; 
F(m+1) = 1;                     % Makes maximum value exactly one
if U(n) >= 1                    % Prevents improper values of probability U
  U(n) = 1;
end
if U(1) <= 0
  U(1) = 1e-9;
end
f  = rowcopy(F,n);              % n rows of F
u  = colcopy(U,m);              % m columns of U
t  = X*((f(:,1:m) < u)&(u <= f(:,2:m+1)))';

dquanplot.m

Descripción del Código:

dquanplot.m Gráfica como escalera grafica la función cuantil para una variable aleatoria simple $X$ . La gráfica son los valores de $X$ versus la función de distribución $F_X$ .

Contestar

% DQUANPLOT file dquanplot.m  Plot of quantile function for a simple rv
% Version of 7/6/95
% Uses stairs to plot the inverse of FX
X  = input('Enter VALUES for X  ');
PX = input('Enter PROBABILITIES for X  ');
m  = length(X);
F  = [0 cumsum(PX)];
XP = [X X(m)];
stairs(F,XP)
grid
title('Plot of Quantile Function')
xlabel('u')
ylabel('t = Q(u)')
hold on
plot(F(2:m+1),X,'o')          % Marks values at jumps
hold off

dsample.m

Descripción del Código:

dsample.m Calcula una muestra a partir de una distribución discreta, determina las frecuencias relativas de los valores y se compara con las probabilidades reales. La entrada consiste en matrices de valor y probabilidad para $X$ y el tamaño de la muestra $n$ . Una matriz $U$ es determinada por un generador de números aleatorios, y se utiliza la función m dquant para calcular los valores de muestra correspondientes. Se calculan y muestran diversos datos sobre la muestra.

Contestar

% DSAMPLE file dsample.m  Simulates sample from discrete population
% Version of 12/31/95 (Display revised 3/24/97)
% Relative frequencies vs probabilities for
% sample from discrete population distribution
X  = input('Enter row matrix of VALUES  ');
PX = input('Enter row matrix of PROBABILITIES  ');
n  = input('Sample size n  ');
U  = rand(1,n);
T  = dquant(X,PX,U);
[x,fr] = csort(T,ones(1,length(T)));
disp('    Value      Prob    Rel freq')
disp([x; PX; fr/n]')
ex = sum(T)/n;
EX = dot(X,PX);
vx = sum(T.^2)/n - ex^2;
VX = dot(X.^2,PX) - EX^2;
disp(['Sample average ex = ',num2str(ex),])
disp(['Population mean E[X] = ',num2str(EX),])
disp(['Sample variance vx = ',num2str(vx),])
disp(['Population variance Var[X] = ',num2str(VX),])

quanplot.m

Descripción del Código:

quanplot.m Traza la función quantile para una función de distribución $F_X$ . Asume que se ha ejecutado el procedimiento dfsetup o acsetup. Se determina un conjunto adecuado $U$ de valores de probabilidad y se utiliza la función m dquant para determinar los valores correspondientes de la función cuantil. Se trazan los resultados.

Contestar

% QUANPLOT file quanplot.m  Plots quantile function for dbn function
% Version of 2/2/96
% Assumes dfsetup or acsetup has been run
% Uses m-function dquant
X  = input('Enter row matrix of values  ');
PX = input('Enter row matrix of probabilities  ');
h  = input('Probability increment h  ');
U  = h:h:1;
T  = dquant(X,PX,U);
U  = [0 U 1];
Te = X(m) + abs(X(m))/20;
T  = [X(1) T Te];
plot(U,T)             % Plot rather than stairs for general case
grid
title('Plot of Quantile Function')
xlabel('u')
ylabel('t = Q(u)')

qsample.m

Descripción del Código:

qsample.m Simula una muestra para una densidad de población dada. Determina los parámetros de muestra y los parámetros de población aproximados Asume que dfsetup o acsetup ha sido ejecutada. Toma como entrada las matrices de distribución $X, PX$ y el tamaño de la muestra $n$ . Utiliza un generador de números aleatorios para obtener la matriz de probabilidad $U$ y utiliza la función m dquant para determinar la muestra. Asume que dfsetup o acsetup ha sido ejecutada.

Contestar

% QSAMPLE file qsample.m  Simulates sample for given population density
% Version of 1/31/96
% Determines sample parameters 
% and approximate population parameters.
% Assumes dfsetup or acsetup has been run
X  = input('Enter row matrix of VALUES  ');
PX = input('Enter row matrix of PROBABILITIES  ');
n  = input('Sample size n =  ');
m  = length(X);
U  = rand(1,n);
T  = dquant(X,PX,U);
ex = sum(T)/n;
EX = dot(X,PX);
vx = sum(T.^2)/n - ex^2;
VX = dot(X.^2,PX) - EX^2;
disp('The sample is in column vector T')
disp(['Sample average ex = ', num2str(ex),])
disp(['Approximate population mean E(X) = ',num2str(EX),])
disp(['Sample variance vx = ',num2str(vx),])
disp(['Approximate population variance V(X) = ',num2str(VX),])

targetset.m

Descripción del Código:

targetset.m Configuración para la llegada a un conjunto de valores objetivo. Se utiliza junto con el procedimiento m targetrun para determinar el número de ensayos necesarios para visitar $k$ un conjunto específico de valores objetivo. La entrada consiste en las matrices de distribución $X, PX$ y el conjunto especificado $E$ de valores objetivo.

Responder

% TARGETSET file targetset.m  Setup for sample arrival at target set
% Version of 6/24/95
X  = input('Enter population VALUES  ');
PX = input('Enter population PROBABILITIES  ');
ms = length(X);
x = 1:ms;                   % Value indices
disp('The set of population values is')
disp(X);
E  = input('Enter the set of target values  ');
ne = length(E);
e  = zeros(1,ne);
for i = 1:ne
  e(i) = dot(E(i) == X,x);  % Target value indices
end
F  = [0 cumsum(PX)];
A  = F(1:ms);
B  = F(2:ms+1);
disp('Call for targetrun')

targetrun.m

Descripción del Código:

targetrun.m Asume que el target m-file ha proporcionado los datos básicos. La entrada consiste en el número $r$ de repeticiones y el número $k$ de los estados objetivo a visitar. Calcula y muestra diversos resultados.

Responder

% TARGETRUN file targetrun.m Number of trials to visit k target values
% Version of 6/24/95  Rev for Version 5.1 1/30/98
% Assumes the procedure targetset has been run.
r  = input('Enter the number of repetitions  ');
disp('The target set is')
disp(E)
ks = input('Enter the number of target values to visit  ');
if isempty(ks)
  ks = ne;
end
if ks > ne
  ks = ne;
end
clear T             % Trajectory in value indices (reset)
R0 = zeros(1,ms);   % Indicator for target value indices
R0(e) = ones(1,ne);
S  = zeros(1,r);    % Number of trials for each run (reset)
for k = 1:r
  R = R0;
  i = 1;
  while sum(R) > ne - ks
    u = rand(1,1);
    s = ((A < u)&(u <= B))*x';
    if R(s) == 1     % Deletes indices as values reached
      R(s) = 0;
    end
    T(i) = s;
    i = i+1;
  end
  S(k) = i-1;
end
if r == 1
  disp(['The number of trials to completion is ',int2str(i-1),])
  disp(['The initial value is ',num2str(X(T(1))),])
  disp(['The terminal value is ',num2str(X(T(i-1))),])
  N  = 1:i-1;
  TR = [N;X(T)]';
  disp('To view the trajectory, call for TR')
else
  [t,f] = csort(S,ones(1,r));
  D  = [t;f]';
  p  = f/r;
  AV = dot(t,p);
  SD = sqrt(dot(t.^2,p) - AV^2);
  MN = min(t);
  MX = max(t);
  disp(['The average completion time is ',num2str(AV),])
  disp(['The standard deviation is ',num2str(SD),])
  disp(['The minimum completion time is ',int2str(MN),])
  disp(['The maximum completion time is ',int2str(MX),])
  disp(' ')
  disp('To view a detailed count, call for D.')
  disp('The first column shows the various completion times;')
  disp('the second column shows the numbers of trials yielding those times')
  plot(t,cumsum(p))
  grid
  title('Fraction of Runs t Steps or Less')
  ylabel('Fraction of runs')
  xlabel('t = number of steps to complete run')
end

Demanda compuesta

El siguiente patrón proporciona un modelo útil en muchas situaciones. Considerar

$D = \sum_{k = 0}^{N} Y_k$

donde $Y_0 = 0$ , y la clase $\{Y_k: 1 \le k\}$ es iid, independiente de la variable aleatoria de conteo $N$ . Una interpretación natural es considerar $N$ que es el número de clientes en una tienda y $Y_k$ la cantidad comprada por el cliente $k$ th. Entonces $D$ es la demanda total de los clientes reales. De ahí que llamemos a $D$ la demanda compuesta.

gend.m

Descripción del Código:

gend.m Utiliza coeficientes de las funciones generadoras para $N$ y $Y$ para calcular, en el caso entero, la distribución marginal para la demanda compuesta $D$ y la distribución conjunta para $\{N, D\}$

Responder

% GEND file gend.m   Marginal and joint dbn for integer compound demand
% Version of 5/21/97
% Calculates marginal distribution for compound demand D
% and joint distribution for {N,D} in the integer case
% Do not forget zero coefficients for missing powers
% in the generating functions for N, Y
disp('Do not forget zero coefficients for missing powers')
gn = input('Enter gen fn COEFFICIENTS for gN  ');
gy = input('Enter gen fn COEFFICIENTS for gY  ');
n  = length(gn) - 1;           % Highest power in gN
m  = length(gy) - 1;           % Highest power in gY
P  = zeros(n + 1,n*m + 1);     % Base for generating P
y  = 1;                        % Initialization
P(1,1) = gn(1);                % First row of P (P(N=0) in the first position)
for i = 1:n                    % Row by row determination of P
   y  = conv(y,gy);            % Successive powers of gy
   P(i+1,1:i*m+1) = y*gn(i+1); % Successive rows of P
end
PD = sum(P);                   % Probability for each possible value of D
a  = find(gn);                 % Location of nonzero N probabilities
b  = find(PD);                 % Location of nonzero D probabilities
P  = P(a,b);                   % Removal of zero rows and columns
P  = rot90(P);                 % Orientation as on the plane
N  = 0:n;
N  = N(a);                     % N values with positive probabilites
PN = gn(a);                    % Positive N probabilities
Y  = 0:m;                      % All possible values of Y
Y  = Y(find(gy));              % Y values with positive probabilities
PY = gy(find(gy));             % Positive Y proabilities
D  = 0:n*m;                    % All possible values of D
PD = PD(b);                    % Positive D probabilities
D  = D(b);                     % D values with positive probabilities
gD = [D; PD]';                 % Display combination
disp('Results are in N, PN, Y, PY, D, PD, P')
disp('May use jcalc or jcalcf on N, D, P')
disp('To view distribution for D, call for gD')

gendf.m

Descripción del Código:

gendf.m function [d, pd] = gendf (gn, gy) es una versión de función de gend, que permite nombrar arbitrariamente las variables. Calcula la distribución para $D$ , pero no la distribución conjunta para $\{N, D\}$

Responder

function [d,pd] = gendf(gn,gy)
% GENDF [d,pd] = gendf(gN,gY) Function version of gend.m
% Calculates marginal for D in the integer case 
% Version of 5/21/97
% Do not forget zero coefficients for missing powers
% in the generating functions for N, Y
n  = length(gn) - 1;           % Highest power in gN
m  = length(gy) - 1;           % Highest power in gY
P  = zeros(n + 1,n*m + 1);     % Base for generating P
y  = 1;                        % Initialization
P(1,1) = gn(1);                % First row of P (P(N=0) in the first position)
for i = 1:n                    % Row by row determination of P
   y  = conv(y,gy);            % Successive powers of gy
   P(i+1,1:i*m+1) = y*gn(i+1); % Successive rows of P
end
PD = sum(P);                   % Probability for each possible value of D
D  = 0:n*m;                    % All possible values of D
b  = find(PD);                 % Location of nonzero D probabilities
d  = D(b);                     % D values with positive probabilities
pd = PD(b);                    % Positive D probabilities

mgd.m

Descripción del Código:

mgd.m Utiliza coeficientes para la función generadora para $N$ y la distribución para simple $Y$ para calcular la distribución para la demanda compuesta.

Responder

% MGD file mgd.m  Moment generating function for compound demand
% Version of 5/19/97
% Uses m-functions csort, mgsum
disp('Do not forget zeros coefficients for missing')
disp('powers in the generating function for N')
disp(' ') 
g  = input('Enter COEFFICIENTS for gN  ');
y  = input('Enter VALUES for Y  ');
p  = input('Enter PROBABILITIES for Y  ');
n  = length(g);               % Initialization
a  = 0;
b  = 1;
D  = a;
PD = g(1);
for i = 2:n
  [a,b] = mgsum(y,a,p,b);
  D  = [D a];
  PD = [PD b*g(i)];
  [D,PD] = csort(D,PD);
end
r = find(PD>1e-13); 
D = D(r);                     % Values with positive probability
PD = PD(r);                   % Corresponding probabilities
mD = [D; PD]';                % Display details
disp('Values are in row matrix D; probabilities are in PD.')
disp('To view the distribution, call for mD.')

Ejercicio $\PageIndex{1}$

Descripción del Código:

mgdf.m function [d, pd] = mgdf (pn, y, py) es una versión de función de mgd, que permite nombrar arbitrariamente las variables. La matriz de entrada $pn$ es la matriz de coeficientes para la función de generación de variables aleatorias de conteo. Deben incluirse ceros para los poderes faltantes. Las matrices $y, py$ son los valores reales y probabilidades de la variable aleatoria de demanda.

Responder

function [d,pd] = mgdf(pn,y,py)
% MGDF [d,pd] = mgdf(pn,y,py)  Function version of mgD
% Version of 5/19/97
% Uses m-functions mgsum and csort
% Do not forget zeros coefficients for missing
% powers in the generating function for N 
n  = length(pn);               % Initialization
a  = 0;
b  = 1;
d  = a;
pd = pn(1);
for i = 2:n
  [a,b] = mgsum(y,a,py,b);
  d  = [d a];
  pd = [pd b*pn(i)];
  [d,pd] = csort(d,pd);
end
a  = find(pd>1e-13);          % Location of positive probabilities
pd = pd(a);                   % Positive probabilities
d  = d(a);                    % D values with positive probability

Ejercicio $\PageIndex{1}$

Descripción del Código:

randbern.m Que S sea el número de éxitos en un número aleatorio $N$ de ensayos de Bernoulli, con probabilidad $p$ de éxito en cada ensayo. El procedimiento randbern toma como entradas la probabilidad $p$ de éxito y las matrices de distribución $N$ , $PN$ para la variable aleatoria de conteo $N$ y calcula la distribución conjunta para $\{N, S\}$ y la distribución marginal para $S$ .

Responder

% RANDBERN file randbern.m   Random number of Bernoulli trials
% Version of 12/19/96; notation modified 5/20/97
% Joint and marginal distributions for a random number of Bernoulli trials
% N is the number of trials
% S is the number of successes
p  = input('Enter the probability of success  ');
N  = input('Enter VALUES of N  ');
PN = input('Enter PROBABILITIES for N  ');
n  = length(N);
m  = max(N);
S  = 0:m;
P  = zeros(n,m+1);
for i = 1:n
  P(i,1:N(i)+1) = PN(i)*ibinom(N(i),p,0:N(i));
end
PS = sum(P);
P  = rot90(P);
disp('Joint distribution N, S, P, and marginal PS')

Simulación de sistemas Markov

Ejercicio $\PageIndex{1}$

Descripción del Código:

inventory1.m Calcula la matriz de transición para una política de $(m, M)$ inventario. Al final de cada periodo, si el stock es menor que un punto de reorden $m$ , el stock se repone al nivel $M$ . La demanda en cada periodo es una variable aleatoria de valor entero $Y$ . La entrada consiste en los parámetros $m, M$ y la distribución para Y como una simple variable aleatoria (o una aproximación discreta).

Responder

% INVENTORY1 file inventory1.m  Generates P for (m,M) inventory policy
% Version of 1/27/97
% Data for transition probability calculations
% for (m,M) inventory policy
M  = input('Enter value M of maximum stock  ');
m  = input('Enter value m of reorder point  ');
Y  = input('Enter row vector of demand values  ');
PY = input('Enter demand probabilities  ');
states = 0:M;
ms = length(states);
my = length(Y);
% Calculations for determining P
[y,s] = meshgrid(Y,states);
T  =  max(0,M-y).*(s < m) + max(0,s-y).*(s >= m);
P  = zeros(ms,ms);
for i = 1:ms
   [a,b] = meshgrid(T(i,:),states);
   P(i,:) = PY*(a==b)';
end
disp('Result is in matrix P')

ramap.m

Descripción del Código:

branchp.m Calcula la matriz de transición para un proceso de ramificación simple con una población máxima especificada. La entrada consiste en el valor máximo de población $M$ y la matriz de coeficientes para la función generadora para las variables aleatorias de propagación individuales $Z_i$ . Esta última matriz debe incluir coeficientes cero para las potencias faltantes.

Responder

% BRANCHP file branchp.m  Transition P for simple branching process
% Version of 7/25/95
% Calculates transition matrix for a simple branching 
% process with specified maximum population.
disp('Do not forget zero probabilities for missing values of Z')
PZ = input('Enter PROBABILITIES for individuals  ');
M  = input('Enter maximum allowable population  ');
mz = length(PZ) - 1;
EZ = dot(0:mz,PZ);
disp(['The average individual propagation is ',num2str(EZ),])
P  = zeros(M+1,M+1);
Z  = zeros(M,M*mz+1);
k  = 0:M*mz;
a  = min(M,k);
z  = 1;
P(1,1) = 1;
for i = 1:M                 % Operation similar to gend
  z = conv(PZ,z);
  Z(i,1:i*mz+1) = z;
  [t,p] = csort(a,Z(i,:));
  P(i+1,:) = p;
end  
disp('The transition matrix is P')
disp('To study the evolution of the process, call for branchdbn')

juego de cadenas.m

Descripción del Código:

chainset.m Se configura para simulación de cadenas de Markov. Las entradas son la matriz de transición P el conjunto de estados, y un conjunto opcional de estados objetivo. Los procedimientos de generación de cadena enumerados a continuación asumen que este procedimiento se ha ejecutado.

Responder

% CHAINSET file chainset.m Setup for simulating Markov chains
% Version of 1/2/96 Revise 7/31/97 for version 4.2 and 5.1
P  = input('Enter the transition matrix  ');
ms = length(P(1,:));
MS = 1:ms;
states = input('Enter the states if not 1:ms  ');
if isempty(states)
  states = MS;
end
disp('States are')
disp([MS;states]')
PI = input('Enter the long-run probabilities  ');
F  = [zeros(1,ms); cumsum(P')]';
A  = F(:,MS);
B  = F(:,MS+1);
e  = input('Enter the set of target states  ');
ne = length(e);
E  = zeros(1,ne);
for i = 1:ne
  E(i) = MS(e(i)==states);
end
disp(' ')
disp('Call for for appropriate chain generating procedure')

mchain.m

Descripción del Código:

mchain.m Asume que se ha corrido el juego de platos. Genera trayectoria de longitud especificada, con estado inicial especificado.

Responder

% MCHAIN file mchain.m  Simulation of Markov chains
% Version of 1/2/96  Revised 7/31/97 for version 4.2 and 5.1
% Assumes the procedure chainset has been run
n  = input('Enter the number n of stages   ');
st = input('Enter the initial state  ');
if ~isempty(st)
  s  = MS(st==states);
else
  s = 1;
end
T  = zeros(1,n);           % Trajectory in state numbers
U  = rand(1,n);
for i = 1:n
  T(i) = s;
  s = ((A(s,:) < U(i))&(U(i) <= B(s,:)))*MS';
end
N  = 0:n-1;
tr = [N;states(T)]';
n10 = min(n,11);
TR = tr(1:n10,:);
f  = ones(1,n)/n;
[sn,p] = csort(T,f);
if isempty(PI)
  disp('     State     Frac')
  disp([states; p]')
else
  disp('     State     Frac       PI')
  disp([states; p; PI]')
end
disp('To view the first part of the trajectory of states, call for TR')

llegada.m

Descripción del Código:

arrival.m Asume que se ha corrido el juego de platos. Calcula repetidamente la hora de llegada a un conjunto prescrito de estados.

Responder

% ARRIVAL file arrival.m  Arrival time to a set of states
% Version of 1/2/96  Revised 7/31/97 for version 4.2 and 5.1
% Calculates repeatedly the arrival
% time to a prescribed set of states.
% Assumes the procedure chainset has been run.
r  = input('Enter the number of repetitions  ');
disp('The target state set is:')
disp(e)
st = input('Enter the initial state  ');
if ~isempty(st)
  s1 = MS(st==states); % Initial state number
else
  s1 = 1;
end
clear T                % Trajectory in state numbers (reset)
S  = zeros(1,r);       % Arrival time for each rep  (reset)
TS = zeros(1,r);       % Terminal state number for each rep (reset)
for k = 1:r
  R  = zeros(1,ms);    % Indicator for target state numbers
  R(E) = ones(1,ne);   % reset for target state numbers
  s  = s1;
  T(1) = s;
  i  = 1;
  while R(s) ~= 1      % While s is not a target state number
    u = rand(1,1);
    s = ((A(s,:) < u)&(u <= B(s,:)))*MS';
    i = i+1;
    T(i) = s;
  end
  S(k) = i-1;          % i is the number of stages; i-1 is time
  TS(k) = T(i);
end
[ts,ft] = csort(TS,ones(1,r));  % ts = terminal state numbers  ft = frequencies
fts = ft/r;                     % Relative frequency of each ts
[a,at]  = csort(TS,S);          % at = arrival time for each ts
w  = at./ft;                    % Average arrival time for each ts
RES = [states(ts); fts; w]';
disp(' ')
if r == 1
  disp(['The arrival time is ',int2str(i-1),])
  disp(['The state reached is ',num2str(states(ts)),])
  N = 0:i-1;
  TR = [N;states(T)]';
  disp('To view the trajectory of states, call for TR')
else
  disp(['The result of ',int2str(r),' repetitions is:'])
  disp('Term state  Rel Freq   Av time')
  disp(RES)
  disp(' ')
  [t,f]  = csort(S,ones(1,r));  % t = arrival times   f = frequencies
  p  = f/r;                     % Relative frequency of each t
  dbn = [t; p]';
  AV = dot(t,p);
  SD = sqrt(dot(t.^2,p) - AV^2);
  MN = min(t);
  MX = max(t);
  disp(['The average arrival time is ',num2str(AV),])
  disp(['The standard deviation is ',num2str(SD),])
  disp(['The minimum arrival time is ',int2str(MN),])
  disp(['The maximum arrival time is ',int2str(MX),])
  disp('To view the distribution of arrival times, call for dbn')
  disp('To plot the arrival time distribution, call for plotdbn')
end

recurrence.m

Descripción del Código:

recurrence.m Asume que se ha corrido el juego de cadenas. Calcula repetidamente el tiempo de recurrencia a un conjunto prescrito de estados, si el estado inicial está en el conjunto; de lo contrario calcula la hora de llegada.

Responder

% RECURRENCE file recurrence.m Recurrence time to a set of states
% Version of 1/2/96  Revised 7/31/97 for version 4.2 and 5.1
% Calculates repeatedly the recurrence time
% to a prescribed set of states, if initial
% state is in the set; otherwise arrival time.
% Assumes the procedure chainset has been run.
r  = input('Enter the number of repititions  ');
disp('The target state set is:')
disp(e)
st = input('Enter the initial state  ');
if ~isempty(st)
  s1 = MS(st==states); % Initial state number
else
  s1 = 1;
end
clear T                % Trajectory in state numbers (reset)
S  = zeros(1,r);       % Recurrence time for each rep  (reset)
TS = zeros(1,r);       % Terminal state number for each rep (reset)
for k = 1:r
  R  = zeros(1,ms);    % Indicator for target state numbers
  R(E) = ones(1,ne);   % reset for target state numbers
  s  = s1;
  T(1) = s;
  i  = 1;
  if R(s) == 1
    u = rand(1,1);
    s = ((A(s,:) < u)&(u <= B(s,:)))*MS';
    i = i+1;
    T(i) = s;
  end
  while R(s) ~= 1      % While s is not a target state number
    u = rand(1,1);
    s = ((A(s,:) < u)&(u <= B(s,:)))*MS';
    i = i+1;
    T(i) = s;
  end
  S(k) = i-1;          % i is the number of stages; i-1 is time
  TS(k) = T(i);
end
[ts,ft]  = csort(TS,ones(1,r)); % ts = terminal state numbers  ft = frequencies
fts = ft/r;                % Relative frequency of each ts
[a,tt]  = csort(TS,S);    % tt = total time for each ts
w  = tt./ft;               % Average time for each ts
RES = [states(ts); fts; w]';
disp(' ')
if r == 1
  disp(['The recurrence time is ',int2str(i-1),])
  disp(['The state reached is ',num2str(states(ts)),])
  N = 0:i-1;
  TR = [N;states(T)]';
  disp('To view the trajectory of state numbers, call for TR')
else
  disp(['The result of ',int2str(r),' repetitions is:'])
  disp('Term state  Rel Freq   Av time')
  disp(RES)
  disp(' ')
  [t,f]  = csort(S,ones(1,r));  % t = recurrence times   f = frequencies
  p  = f/r;                      % Relative frequency of each t
  dbn = [t; p]';
  AV = dot(t,p);
  SD = sqrt(dot(t.^2,p) - AV^2);
  MN = min(t);
  MX = max(t);
  disp(['The average recurrence time is ',num2str(AV),])
  disp(['The standard deviation is ',num2str(SD),])
  disp(['The minimum recurrence time is ',int2str(MN),])
  disp(['The maximum recurrence time is ',int2str(MX),])
  disp('To view the distribution of recurrence times, call for dbn')
  disp('To plot the recurrence time distribution, call for plotdbn')
end

kvis.m

Descripción del Código:

kvis.m Asume que se ha corrido el juego de platos. Calcula repetidamente el tiempo para completar las visitas a un determinado $k$ de los estados en un conjunto prescrito.

Responder

% KVIS file kvis.m  Time to complete k visits to a set of states
% Version of 1/2/96 Revised 7/31/97 for version 4.2 and 5.1
% Calculates repeatedly the time to complete
% visits to k of the states in a prescribed set.
% Default is visit to all the target states.
% Assumes the procedure chainset has been run.
r  = input('Enter the number of repetitions  ');
disp('The target state set is:')
disp(e)
ks = input('Enter the number of target states to visit  ');
if isempty(ks)
  ks = ne;
end
if ks > ne
  ks = ne;
end
st = input('Enter the initial state  ');
if ~isempty(st)
  s1 = MS(st==states); % Initial state number
else
  s1 = 1;
end
disp(' ')
clear T                    % Trajectory in state numbers (reset)
R0 = zeros(1,ms);          % Indicator for target state numbers
R0(E) = ones(1,ne);        % reset
S = zeros(1,r);            % Terminal transitions for each rep (reset)
for k = 1:r
  R = R0;
  s  = s1;
  if R(s) == 1
    R(s) = 0;
  end
  i  = 1;
  T(1) = s;
  while sum(R) > ne - ks
    u = rand(1,1);
    s = ((A(s,:) < u)&(u <= B(s,:)))*MS';
    if R(s) == 1
      R(s) = 0;
    end
    i = i+1;
    T(i) = s;
  end
  S(k) = i-1;
end
if r == 1
  disp(['The time for completion is ',int2str(i-1),])
  N = 0:i-1;
  TR = [N;states(T)]';
  disp('To view the trajectory of states, call for TR')
else
  [t,f]  = csort(S,ones(1,r));
  p  = f/r;
  D  = [t;f]';
  AV = dot(t,p);
  SD = sqrt(dot(t.^2,p) - AV^2);
  MN = min(t);
  MX = max(t);
  disp(['The average completion time is ',num2str(AV),])
  disp(['The standard deviation is ',num2str(SD),])
  disp(['The minimum completion time is ',int2str(MN),])
  disp(['The maximum completion time is ',int2str(MX),])
  disp(' ')
  disp('To view a detailed count, call for D.')
  disp('The first column shows the various completion times;')
  disp('the second column shows the numbers of trials yielding those times')
end

plotdbn

Descripción del Código:

plotdbn Se utiliza después de la llegada o recurrencia de los m-procedimientos para trazar la distribución del tiempo de llegada o

Responder

% PLOTDBN file plotdbn.m 
% Version of 1/23/98
% Plot arrival or recurrence time dbn
% Use after procedures arrival or recurrence 
% to plot arrival or recurrence time distribution
plot(t,p,'-',t,p,'+')
grid
title('Time Distribution')
xlabel('Time in number of transitions')
ylabel('Relative frequency')

Características de MATLAB

Bloques de construcción auxiliares definidos por el usuario

Vectores y probabilidades minterm

Eventos independientes

Probabilidad condicional e idependencia condicional

Bernoulli y ensayos multinomiales

Algunos problemas de coincidencia

Distribuciones

Dstribuciones binomiales, de Poisson y Gaussianas

Configuración para variables aleatorias simples

Configuración para variables aleatorias generales

Configuración para variables aleatorias simples independientes

Cálculos para variables aleatorias

Cálculos y pruebas para variables aleatorias independientes

Funciones cuantiles para distribuciones delimitadas

Demanda compuesta

Simulación de sistemas Markov

Support Center

How can we help?