risk-assessment

What this skill does

# Risk Assessment Skill

```yaml
name: risk-assessment
version: 1.0.0
category: sme
tags: [risk, probabilistic, monte-carlo, reliability, uncertainty, sensitivity-analysis, marine-safety]
created: 2026-01-06
updated: 2026-01-06
author: Claude
description: |
  Expert risk assessment and probabilistic analysis for marine and offshore
  operations. Includes Monte Carlo simulations, reliability calculations,
  uncertainty quantification, and decision-making under risk.
```

## When to Use This Skill

Use this skill when you need to:
- Perform Monte Carlo simulations for uncertainty quantification
- Calculate system reliability and failure probabilities
- Conduct sensitivity analysis to identify critical parameters
- Create risk matrices for hazard assessment
- Perform probabilistic design and analysis
- Quantify uncertainties in marine operations
- Make decisions under uncertainty with risk metrics
- Validate designs against reliability targets

## Core Knowledge Areas

### 1. Monte Carlo Simulation

Basic Monte Carlo framework:

```python
import numpy as np
from scipy import stats
from dataclasses import dataclass
from typing import Callable, Dict, List, Tuple, Optional
import pandas as pd

@dataclass
class RandomVariable:
    """Statistical distribution for a random variable."""
    name: str
    distribution: str  # 'normal', 'lognormal', 'uniform', 'weibull', etc.
    parameters: dict  # Distribution-specific parameters

    def sample(self, size: int = 1) -> np.ndarray:
        """
        Generate random samples from distribution.

        Args:
            size: Number of samples

        Returns:
            Array of random samples

        Example:
            >>> rv = RandomVariable(
            ...     name='wave_height',
            ...     distribution='weibull',
            ...     parameters={'c': 2.0, 'scale': 3.5}
            ... )
            >>> samples = rv.sample(1000)
        """
        if self.distribution == 'normal':
            return np.random.normal(
                loc=self.parameters['mean'],
                scale=self.parameters['std'],
                size=size
            )
        elif self.distribution == 'lognormal':
            return np.random.lognormal(
                mean=self.parameters['mean'],
                sigma=self.parameters['std'],
                size=size
            )
        elif self.distribution == 'uniform':
            return np.random.uniform(
                low=self.parameters['min'],
                high=self.parameters['max'],
                size=size
            )
        elif self.distribution == 'weibull':
            return stats.weibull_min.rvs(
                c=self.parameters['c'],
                scale=self.parameters['scale'],
                size=size
            )
        elif self.distribution == 'exponential':
            return np.random.exponential(
                scale=self.parameters['scale'],
                size=size
            )
        else:
            raise ValueError(f"Unknown distribution: {self.distribution}")

def monte_carlo_simulation(
    model: Callable,
    random_variables: List[RandomVariable],
    n_samples: int = 10000,
    seed: Optional[int] = None
) -> Dict[str, np.ndarray]:
    """
    Perform Monte Carlo simulation.

    Args:
        model: Function that takes random variable samples and returns result
        random_variables: List of RandomVariable objects
        n_samples: Number of Monte Carlo samples
        seed: Random seed for reproducibility

    Returns:
        Dictionary with input samples and output results

    Example:
        >>> def mooring_tension_model(wave_height, current_speed):
        ...     # Simplified mooring tension model
        ...     return 1000 + 150 * wave_height + 50 * current_speed
        >>>
        >>> rv_wave = RandomVariable(
        ...     'wave_height',
        ...     'weibull',
        ...     {'c': 2.0, 'scale': 3.5}
        ... )
        >>> rv_current = RandomVariable(
        ...     'current_speed',
        ...     'normal',
        ...     {'mean': 1.0, 'std': 0.3}
        ... )
        >>>
        >>> results = monte_carlo_simulation(
        ...     model=lambda samples: mooring_tension_model(
        ...         samples['wave_height'],
        ...         samples['current_speed']
        ...     ),
        ...     random_variables=[rv_wave, rv_current],
        ...     n_samples=10000
        ... )
    """
    if seed is not None:
        np.random.seed(seed)

    # Generate samples for all random variables
    samples = {}
    for rv in random_variables:
        samples[rv.name] = rv.sample(n_samples)

    # Run model for all samples
    outputs = model(samples)

    # Combine inputs and outputs
    results = {**samples, 'output': outputs}

    return results

def calculate_statistics(
    data: np.ndarray,
    percentiles: List[float] = [5, 50, 95]
) -> dict:
    """
    Calculate statistical parameters from Monte Carlo results.

    Args:
        data: Array of simulation results
        percentiles: Percentile values to calculate

    Returns:
        Dictionary with statistics

    Example:
        >>> stats = calculate_statistics(results['output'])
        >>> print(f"Mean: {stats['mean']:.1f}")
        >>> print(f"Std: {stats['std']:.1f}")
        >>> print(f"95th percentile: {stats['p95']:.1f}")
    """
    stats_dict = {
        'mean': np.mean(data),
        'std': np.std(data),
        'min': np.min(data),
        'max': np.max(data),
        'median': np.median(data),
        'cv': np.std(data) / np.mean(data)  # Coefficient of variation
    }

    # Add percentiles
    for p in percentiles:
        stats_dict[f'p{p}'] = np.percentile(data, p)

    return stats_dict
```

### 2. Reliability Analysis

Calculate reliability and failure probability:

```python
def calculate_reliability(
    response_data: np.ndarray,
    limit_state: float,
    mode: str = 'less_than'
) -> dict:
    """
    Calculate reliability from Monte Carlo results.

    Args:
        response_data: Array of response values
        limit_state: Limit state threshold
        mode: 'less_than' or 'greater_than'

    Returns:
        Dictionary with reliability metrics

    Example:
        >>> # Mooring tension should be < 8000 kN
        >>> reliability = calculate_reliability(
        ...     response_data=results['output'],
        ...     limit_state=8000,
        ...     mode='less_than'
        ... )
        >>> print(f"Reliability: {reliability['reliability']:.4f}")
        >>> print(f"Probability of failure: {reliability['pf']:.6f}")
    """
    n_samples = len(response_data)

    if mode == 'less_than':
        # Failure when response >= limit_state
        failures = response_data >= limit_state
    elif mode == 'greater_than':
        # Failure when response <= limit_state
        failures = response_data <= limit_state
    else:
        raise ValueError(f"Unknown mode: {mode}")

    n_failures = np.sum(failures)
    pf = n_failures / n_samples  # Probability of failure
    reliability = 1 - pf

    # Reliability index (beta)
    # β = -Φ^(-1)(Pf) where Φ is standard normal CDF
    if pf > 0 and pf < 1:
        beta = -stats.norm.ppf(pf)
    elif pf == 0:
        beta = np.inf
    else:  # pf == 1
        beta = -np.inf

    return {
        'n_samples': n_samples,
        'n_failures': n_failures,
        'pf': pf,
        'reliability': reliability,
        'beta': beta,
        'target_met': reliability >= 0.9999  # Example: 4-9's reliability
    }

def system_reliability_series(
    component_reliabilities: List[float]
) -> float:
    """
    Calculate system reliability for series system.

    Series system: All components must work for system to work.

    Args:
        component_reliabilities: List of component reliabilities

    Returns:
        System reliability

    Example:
        >>> # 3 mooring lines in series (all must hold)
        >>> R_sys = system_reliability_series([0.999, 0.9995