En esta tesis, se proponen procedimientos Bayesianos para la inferencia, predicción y diseño de diferentes sistemas de colas que incluyen modelos con distribuciones generales para el tiempo entre las llegadas y/o el tiempo de servicio, con uno o varios servidores y con capacidad finita e infinita. Con las observaciones tomadas del proceso de llegadas y de servicio en el sistema, en primer lugar, se desarrollan métodos de estimación Bayesiana de densidades que permiten aproximar las distribuciones desconocidas asociadas a estos procesos. La estimación está basada en mixturas de distribuciones con un número desconocido de componentes que se aborda con diferentes métodos MCMC de dimensión paramétrica variable. Los modelos de mixtura seleccionados son de tipo PH, lo cual ofrece muy buenas propiedades para su posterior aplicación en colas. A continuación, se analiza la congestión del tráfico en cada sistema observado basándose en la inferencia desarrollada sobre los parámetros del mismo. En caso de que el sistema sea estable, se describen procedimientos para la estimación Monte Carlo de las medidas de interés de cada modelo en equilibrio, tales como las distribuciones predictivas del número de clientes en el sistema, del tiempo de espera en cola o de la longitud de los periodos de ocupación. La predicción Bayesiana desarrollada se basa en diferentes propiedades de los sistemas en los que intervienen distribuciones PH y en algunos resultados conocidos de la Teoría de Colas clásica. En los modelos de colas con varios servidores, se abordan también problemas de diseño en los que el objetivo consiste en decidir el número óptimo de servidores que minimizan una función de coste que depende de las distribuciones estacionarias estimadas. La metodología desarrollada se ilustra con datos reales procedentes de un hospital geriátrico de Londres y un establecimiento bancario de Madrid. Por último, se analiza el comportamiento transitorio y el periodo de ocupación de modelos de colas completamente generales con un único servidor. Se obtiene un resultado que permite extraer numéricamente las raíces complejas de unas ecuaciones implicadas en las distribuciones transitorias de interés.